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Abstract Locally adaptive differential frames (gauge 
frames) are a well-known effective tool in image analy¬ 
sis, used in differential invariants and PDE-flows. How¬ 
ever, at complex structures such as crossings or junc¬ 
tions, these frames are not well-defined. Therefore, we 
generalize the notion of gauge frames on images to gauge 
frames on data representations U : xi -A M de¬ 

fined on the extended space of positions and orienta¬ 
tions, which we relate to data on the roto-translation 
group SE{d ), d = 2,3. This allows to define multiple 
frames per position, one per orientation. We compute 
these frames via exponential curve fits in the extended 
data representations in SE(d). These curve fits mini¬ 
mize first or second order variational problems which 
are solved by spectral decomposition of, respectively, a 
structure tensor or Hessian of data on SE(d). We in¬ 
clude these gauge frames in differential invariants and 
crossing preserving PDE-flows acting on extended data 
representation U and we show their advantage com¬ 
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pared to the standard left-invariant frame on SE(d). 
Applications include crossing-preserving filtering and 
improved segmentations of the vascular tree in retinal 
images, and new 3D extensions of coherence-enhancing 
diffusion via invertible orientation scores. 

Keywords Roto-Translation Group • Gauge Frames • 
Exponential Curves • Nonlinear Diffusion • Left- 
invariant Image Processing • Orientation Scores 


1 Introduction 

Many existing image analysis techniques rely on differ¬ 
ential frames that are locally adapted to image data. 
This includes methods based on differential invariants 
[63l[4Hl34ll52] . partial differential equations [6311741140] . 
and non-linear and morphological scale spaces ns 
US, used in various image processing tasks such as track¬ 
ing and line detection [6!, corner detection and edge fo¬ 
cussing mm , segmentation [69] . active contours on 
S3, DTI data processing [481147] , feature based cluster¬ 
ing etc. These local coordinate frames (also known as 
‘gauge frames’ according to Emun]) provide differen¬ 
tial frames directly adapted to the local image structure 
via a structure tensor or a Hessian of the image. Typi¬ 
cally the structure tensor (based on 1st order Gaussian 
derivatives) is used for adapting to edge-like structures 
while the Hessian (based on 2nd order Gaussian deriva¬ 
tives) is used for adapting to line-like structures. The 
primary benefit of the gauge frames is that they allow 
to include adaptation for anisotropy and curvature in 
a rotation and translation invariant way. See Fig. |TJ 
where we have depicted local adaptive frames based on 
eigenvector decomposition of the image Hessian at some 
given scale, of the MR-image in the background. 
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Fig. 1: Left: Locally adaptive frames (gauge frames) 
in the image domain computed as the eigenvectors of 
the Hessian of the image at each location. Right: Such 
gauge frames can be used for adaptive anisotropic dif¬ 
fusion and geometric reasoning. However, at complex 
structures such as blob-structures/crossings, the gauge 
frames are ill-defined causing fluctuations. 


It is sometimes problematic that such locally adapted 
differential frames are directly placed in the image do¬ 
main R d (d = 2,3), as at the vicinity of complex struc¬ 
tures, e.g. crossings, textures, bifurcations, one typi¬ 
cally requires multiple local spatial coordinate frames. 
To this end, one effective alternative is to extend the 
image domain to the joint space of positions and ori¬ 
entations xi 5' d_1 . The advantage is that it allows 
to disentangle oriented structures involved in crossings, 
and to include curvature, cf. Fig.[2j Such extended do¬ 
main techniques rely on various kinds of lifting, such as 
coherent state transforms (also known as invertible ori¬ 
entation scores) P1I23M1I3K] . continuous wavelet trans¬ 
forms [24l[28ll66l l6] . orientation lifts [76lfl2] . or orienta¬ 
tion channel representations |33j. In case one has to 
deal with more complex diffusion weighted MRI tech¬ 
niques, the data in extended position orientation do¬ 
main can be obtained after a modelling procedure as 
in [?nirrnmifi8| . In this article we will not discuss in 
detail on how such a new image representation or lift 
U : R d x S d ~ x R is to be constructed from grey-scale 

image / : R d R, and we assume it to be a sufficiently 
smooth given input. Here V (x, n) is to be considered as 
a probability density of finding a local oriented struc¬ 
ture (i.e. an elongated structure) at position x E 
with orientation n E S d ~ x . 

When processing data in the extended position ori¬ 
entation domain it is often necessary to equip the do¬ 
main with a structure that links the data across dif¬ 
ferent orientation channels, in such a way that a no¬ 
tion of alignment between local orientations is taken 
into account. This is achieved by relating data on po¬ 
sitions and orientations to data on the roto-translation 
group SE(d) = R d xi SO(d). This idea resulted in con- 



Fig. 2: We aim for adaptive anisotropic diffusion of 
images which takes into account curvature. At areas 
with low orientation confidence (in blue) isotropic diffu¬ 
sion is required, whereas at areas with high orientation 
confidence (in red) anisotropic diffusion with curvature 
adaptation is required. Application of locally adaptive 
frames in the image domain suffers from interference 
(3rd column), whereas application of locally adaptive 
frames in the domain x S d ~ 1 allows for adaptation 
along all the elongated structures (4th column). 


textual image analysis methods 
[73ll3Tl[T8ll67| and appears in models of low level visual 
perception and their relation with the functional ar- 

Following the conventions in m we denote functions 
on the coupled space of positions and orientations by 
U :M. d x i S d ~ 1 — R. Then, its extension U : SE(d ) —>• R 
is given by: 


tf(x,R) := U (x, Ra) (1) 

for all x E and all rotations R E SO(d ), and given 
reference axis a E S^ -1 . Throughout this article a is 
chosen as follows: 


d = 2 =x a = (1,0) T , d = 3 =X a = (0,0,1) T . (2) 

Then, we can identify the joint space of positions and 
orientations xi S d ~ x by: 

R d x gd-t SE(d)/{{ 0} x SO(d - 1)), (3) 

where this quotient structure is due to ([]]), and where 
SO(d— 1) is identified with all rotations on that map 
reference axis a onto itself. Note that in Eq. 0 the tilde 
indicates we consider data on the group instead of data 
on the quotient. If d = 2 the tildes can be ignored as 
M 2 x S 1 = SE( 2). However, for d > 3 this distinction 
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is crucial and necessary details on © will follow in the 
beginning of Section [6] 

In this article, our quest is to find locally optimal 
differential frames in SE(d) relying on similar Hessian- 
and/or structure-tensor type of techniques for gauge 
frames on images, recall Fig. [T] Then, the frames can 
be used to construct crossing-preserving differential in¬ 
variants and adaptive diffusions of data in SE(d). In 
order to find these optimal frames our main tool is the 
theory of curve fits. Early works on curve fits have been 
presented in [57 where the notion of curvature consis¬ 
tency is applied to inferring local curve orientations, 
based on neighbourhood co-circularity continuation cri¬ 
teria. This approach was extended to 2D texture flow 
inference in [8], by lifting images in position and ori¬ 
entation domain and inferring multiple Cartan frames 
at each point. Our work is embedded in a Lie group 
framework where we consider the notion of exponential 
curve fits via formal variational methods. Exponential 
curves in the SE(d)- curved geometry are the equiva¬ 
lents of straight lines in the Euclidean geometry. If 
d = 2, the spatial projection of these exponential curves 
are osculating circles, which are used: for constructing 
the curvature consistency in [S3, for defining the tensor 
voting fields in [54], and for local modeling association 
fields in [20] . If d = 3, the spatial projection of expo¬ 
nential curves are spirals with constant curvature and 
torsion. Based on co-helicity principles, similar spirals 
have been used in neuroimaging applications [64 or for 
modelling heart fibers [65]. In these works curve fits are 
obtained via efficient discrete optimization techniques, 
which are beyond the scope of this article. 

In Fig. [3] we present an example for d = 2 of the 
overall pipeline of including locally adaptive frames in 
a suitable diffusion operators acting in the lifted do¬ 
main M 2 x S 1 . For d > 2 the same pipeline applies. 
Here, an exponential curve fit 7^ (t) (in blue, with spa¬ 
tial projection in red) at a group element g G SE(d) is 
characterized by (g,c*(g)), i.e. a starting point g and 
an tangent vector c*(<7) that should be aligned with the 
structures of interest. In essence, this paper explains in 
detail how to compute c*(<7) as this will be the princi¬ 
pal direction the differential frame will be aligned with, 
and then gives appropriate conditions for fixing the re¬ 
maining directions in the frame. 

The main contribution of this article is to provide a 
general theory for finding locally adaptive frames in the 
roto-translation group SE(d ), for d = 2,3. Some pre¬ 
liminary work on exponential curve fits of the second 
order on SE( 2) has been presented in [35ll36ll66] . In this 
paper we formalize these previous methods (Theorems 


1 Exponential curves are auto-parallels w.r.t. ‘-’Cartan con¬ 

nection, see Appendix A, Eq. (132). 


[2] and [3]) and we extend them to first-order exponen¬ 
tial curve fits (Theorem [l]). Furthermore, we generalize 
both approaches to the case d = 3 (Theorems |4j [5j 
[6j [7] and [8]). All theorems contain new results except 
for Theorems [2] and [3] The key ingredient is to con¬ 
sider the fits as formal variational curve optimization 
problems with exact solutions derived by spectral de¬ 
composition of structure tensors and Hessians of the 
data U on SE(d). In the SE( 3)-case we show that in 
order to obtain torsion-free exponential curve fits with 
well-posed projection on M 3 xi S' 2 , one must resign to a 
two-fold optimization algorithm. To show the potential 
of considering these locally adaptive frames, we em¬ 
ploy them in medical image analysis applications, in 
improved differential invariants and improved crossing¬ 
preserving diffusions. Here, we provide for the first time 
coherence enhancing diffusions via 3D invertible orien¬ 
tation scores mui, extending previous methods [35l 
1361166] to the 3D Euclidean motion group. 


1.1 Structure of the Article 

We start the body of this article reviewing preliminary 
differential geometry tools in Section [2] Then, in Sec¬ 
tion [3] we describe how a given exponential curve fit 
induces the locally adaptive frame. In Section [4] we 
provide an introduction by reformulating the standard 
gauge frames construction in images in a group theo¬ 
retical setting. This gives a roadmap towards SE( 2)- 
extensions explained in Section [5j where we deal with 
exponential curve fits of the 1st order in Subsection |5.2| 
computed via a structure tensor, and exponential curves 
fits of 2nd order in Section |5.3| computed via the Hes¬ 
sian of the data U. In the latter case we have 2 options 
for the curve optimization problem, one solved by the 
symmetric sum, and one by the symmetric product of 
the non-symmetric Hessian. The curve fits in SE( 2) in 
Section [3J are extended to curve fits in SE( 3) in Sec¬ 
tion [3] It starts with preliminaries on the quotient © 
and then it follows the same structure as the previ¬ 
ous section. Here we present the two-fold algorithm for 
computing the torsion free exponential curve fits. 

In Section[7]we consider experiments regarding med¬ 
ical imaging applications and feasibility studies. We 
first recall the theory of invertible orientation scores 
needed for the applications. In the SE( 2)-case we present 
crossing-preserving multi-scale vessel enhancing filters 
in retinal imaging, and in the 5'F(3)-case we include a 
proof of concept of crossing-preserving (coherence en¬ 
hancing diffusion) steered by gauge frames via invert¬ 
ible 3D orientation scores. 

Finally, there are 5 appendices. Appendix [A] supple¬ 
ments Section [3] by explaining the construction of the 






4 


Duits, Janssen, Hannink, Sanguinetti 


Image / Score U 



Fig. 3: The overall pipeline of image processing f Yf via left-invariant operators In this pipeline we construct 
an invertible orientation score W^f (Section we fit an exponential curve (Section |5|6| ), we obtain the gauge 
frame (Section [ 3 ] and App. |A|), we construct a non-linear diffusion, and finally we apply reconstruction (Section 


7.1). The main focus of this paper is curve fitting, where we compute per element g = {x,y,0) an exponential 
curve fit 7 ® (t) (in blue, with spatial projection in red) with tangent 7 ^ (0) = c *(g) = (c 1 ,c 2 ,c 3 ) T at g. Based on 


this fit we construct for each g a local frame {Bi\ g . 
^ is a non-linear diffusion operator). 


S 2 I 0 , Bs\ q } which are used in our operators on the lift (here 
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frame for d = 2,3. Appendix |B| describes the geometry 
of neighboring exponential curves needed for formulat¬ 
ing the variational problems. Appendix [C] complements 
the two-fold approach in Section [6j Appendix [D] pro¬ 
vides the definition of the Hessian used in the paper. 
Finally, Appendix 1 contains a list of symbols, their 
explanations and references to the equation in which 
they are defined. We advise the reader to keep track 
of this table. Especially, in the more technical sections: 
Section [5] and [6j 


2 Differential Geometrical Tools 

Relating our data to data on the Euclidean motion 
group, via Eq. 0, allows us to use tools from Lie group 
theory and differential geometry. In this section we ex¬ 
plain these tools that are important for our notion of an 
exponential curve fit to smooth data U : SE(d ) M. 
Often, we consider the case d = 2 for basic illustration. 
Later on, in Section [6j we consider the case d = 3 and 
extra technicalities on the quotient structure will enter. 


2.1 The Roto-Translation Group 

The data U : SE(d) M is defined on the group SE(d) 
of rotations and translations acting on M d . As the con¬ 
catenation of two rigid body motions is again a rigid 
body motion, the group SE(d) is equipped with the 
following group product: 

99 ' = (x, R)(x', R') = (Rx' + X, RR'), (d) 

with g = (x, R), g r = (x^R 7 ) G SE(d ), ^ 

where we recognize the semi-direct product structure 
SE(d) = x SO(d ), of the translation group R d with 
rotation group SO(d) = {R G M dxd | R T = R _ * det R= 1}. 
The groups SE(d) and SO(d) have dimension 

r d :=dim(SO(d)) = ^^, 

n d := dim (SE(d)) = = d + r d . 

Note that n 2 = 3, 713 = 6 . One may represent elements 
g from SE(d) by the following matrix representation 

g = M(g) = ^qT 5 which indeed satisfies ^ 
M(gg') = M(g)M(g'). 
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We will often avoid this embedding into the set of in¬ 
vertible (d + 1 ) x (d + 1 ) matrices, in order to focus on 
the geometry rather than the algebra. 

2.2 Left-Invariant Operators 


A ^ 

where A i-a e A = ^ denotes the matrix exponen- 
k=0 

tial from Lie algebra T t (SE(d)) to Lie group SE(d). 
The differential operators {A}™=;l induce a correspond¬ 
ing dual frame which is a basis for the co¬ 

tangent bundle T*(SE(d)). This dual frame is given by 


In image analysis applications operators U T?{U) 
need to be left-invariant and not right-invariant [25, 
135] . Left-invariant operators d? in the extended domain 
correspond to rotation and translation invariant opera¬ 
tors T in the image domain, which is a desirable prop¬ 
erty. On the other hand, right-invariance boils down to 
isotropic operators in the image domain which is an 
undesirable restriction. By definition d> is left-invariant 
and not right-invariant if it commutes with the left- 
regular representation C (and not with the right-regular 
representation TV). Representations C,7Z are given by 

(ChU)(g) = U(h~ 1 g), ( K h U)(g ) = U(gh), (7) 

for all h^g G SE(d). So operator d> must satisfy d>oC g = 
C g od> and 3 o !Z g ^ 1Z g o 3 for all g G SE(d). 


2.3 Left-Invariant Vector Fields and Dual Frame 


A special case of left-invariant operators are left-invariant 
derivatives. More precisely (see Remark [l] below), we 
need to consider left-invariant vector fields g A gi 
as the left-invariant derivative A g depends on the loca¬ 
tion g where it is attached. Intuitively, the left-invariant 
vector fields {Ai }™! L provide a local moving frame of 
reference in the tangent bundle T(SE(d )), that comes 
in naturally when including alignment of local orienta¬ 
tions in the image processing of U. 

Formally, the left-invariant vector fields are obtained 
by taking a basis {A*}™^ G T t (SE(d)) in the tangent 
space at the unity element e := (0,/), and then one 
uses the push-forward (L g )* of the left multiplication 

L g h = gh , ( 8 ) 

to obtain the corresponding tangent vectors in the tan¬ 
gent space T g (SE(d)). Thus one associates to each Ai 
a left-invariant field Ai given by 


Ai\ g = {L g )*Ai, for all g £ SE(d), i = l,...,n d , (9) 

where we consider each Ai as a differential operator on 
smooth locally defined functions (/) given by 


AL 0 - (Lg) *Airf) :- Ai(4> 0 Lg). 


An explicit way to construct and compute the differen¬ 
tial operators Ai\ from A = A| c is via 


A|A = A <i>(g) = lim 

y e —>-0 


< Kge eAi ) - 4>(g) 


(10) 


= 5) for i,j = l,...n d , ( 11 ) 

where 5 l - denotes the Kronecker delta. Then the deriva¬ 
tive of a differentiable function <j) : SE(d) -A M is ex¬ 
pressed as follows 

n d 

d 4 = '%2Ai$w i eT*{SE(d)). (12) 

i= 1 


Remark 1 In differential geometry, there exist two equiv¬ 
alent viewpoints [3j Ch. 2] on tangent vectors A g G 
T g (SE(d )): either one considers them as tangents to lo¬ 
cally defined curves; or one considers them as differen¬ 
tial operators on locally defined functions. The connec¬ 
tion between these viewpoints is as follows. We identify 
a tangent vector q(t) G T^^(SE(d)) with the differen¬ 
tial operator (q(t))( 0 ) := J^ 0 ( 7 (£)) for all locally de¬ 
fined, differentiable, real-valued functions 0 . 


Next we express tangent vectors explicitly in the left- 
invariant moving frame of reference, by taking a direc¬ 
tional derivative: 


g _ . Ud 

— 0(7(t)) = (d<^(7(f)),7(l)) = J^f(t) Ai\ m 4> 

i= 1 


(13) 


n d . . 

with 7 (t) = y; fit) Ai\ W) , and with </> smooth and 
i= 1 

defined on an open set around j(t). Eq. will play 
a crucial role in Section [5] (exponential curve fits for 
d — 2 ) and Section [ 6 ] (exponential curve fits for d = 3). 


Example 1 For d = 2 we take Ai = d x \ t , A 2 = d y \ , 
A 3 = 0q\ . Then we have the left-invariant vector fields 


M {x ,y,9) ; =COS0 ^\ {xy 9) +smO §- y ^ 

M(x,v,6) ■= ~ sind {x ,y,9)+ COsd ly 
A'i\( Xt y,0) ■— d$\(x,y,0)' 

The dual frame is given by 

lj 1 = cos 0 dx + sin Qdy, 
uj 2 = — sin 6 dx + cos Ody, 
cj 3 = dO. 


(14) 


(15) 


For explicit formulas for left-invariant vector fields in 
SE( 3) see [T91I3T] . 


e 
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2.4 Exponential Curves in SE{d ) 


2.5 Left-Invariant Metric Tensor on SE(d) 


Let (cW.cW ) 7 E M d+rd = M nd be a given column 
vector, where = (c 1 ,..., c d ) E denotes the spa¬ 
tial part and c^ 2 ) = (c d+1 ,... , c nd ) E W' d denotes the 

rotational part. The unique exponential curve passing 

n d 

through# E SE(d ) with initial velocity c(g) = ^2 c 1 Ai\ 

i=X 

equals 



(16) 


We use the following (left-invariant) metric tensor: 

d n d 

^^( 7 , 7 )=m 2 Ei^I 2+ E l^l 2 ’ ( 23 ) 

i=l i=d +1 

where 7 = and with stiffness parameter 

(i along any smooth curve 7 in SE(d). Now, for the 
special case of exponential curves, one has 7 * = c 1 is 
constant. The metric allows us to normalize the speed 
along the curves by imposing a normalization constraint 


with Ai = Ai\ denoting a basis of T t (SE(d)). In fact 
such exponential curves satisfy 

n d 

ftf) = ( 17 ) 

2=1 

and thereby have constant velocity in the moving frame 
of reference, i.e. 7* = c 1 in Eq. ([ 13 J) . A way to compute 
the exponentials is via matrix exponentials and & 

Example 2 If d = 2 we have exponential curves: 

%(t) = go e *(c%+^2+c 3 A 3 ) = ,,(*), (9(4)), ( 18 ) 

which are circular spirals with 

x{t) = x 0 + ^(sin(c 3 t+6>o) - sin(< 9 0 )) 

+ ^(cos(c 3 t+6> 0 ) -cos(6> 0 )) , 

2/W = 2/o — ^(cos(c 3 t+6> 0 ) - cos(< 9 0 )) ( 19 ) 

-h^(sin(c 3 t+6> 0 ) - sin(6> 0 )) , 

6>(t) = 6> 0 + tc 3 , 

for the case c 3 7^ 0, and allt > 0 and straight lines with 


x(t) = xq + ^(c 1 cos^o — c 2 sin ^ 0)5 

y(t) = #0 + ^(c 1 sin # 0 + c 2 cos 0 O )> ( 20 ) 

0 (*)« 0 o, 

for the case c 3 = 0, where #0 = (%o,yo,@o) E SE{2). 
See the left panel in Fig. [4] 


Example 3 For d = 3, the formulae for exponential 
curves in SE( 3) are given in for example [TOllIU] . Their 
spatial part are circular spirals with torsion r(t) = 

lc (1) -c (2) l / X ! 

1 || c (i)|i K(t) and curvature 


«(*) = 
+ 


||c (1) | 
sin(t || < 


(cos(t||c^ 2 ^ ||) c^ 2 ^ X c^ 1 ) 


c(2) II 


ill r (2) 


X c 


( 2 ) 


X c 


( 1 ) 


( 21 ) 


Note that their magnitudes are constant: 


Ilct 1 ) x c( 2 )|| 

llc^ll 2 


and |r| 


| c ( 1 ) • c^ 2 ^ | • \k\ 


( 22 ) 


d 


l|c || 2 := I|M m c|| 2 =# 2 £|u 


n d 


2 - ff 2 Y' \A\2 _L ^ \C 1 \ 2 

= 1 i = d-\- 1 

= M 2 ||c( 1 >|| 2 +||c ( 2 )|| 2 = l, 


with := 


yh 0 

0 I r , 


(24) 


We will use this constraint in the fitting procedure in 
order to ensure that our exponential curves © are 
parameterized by Riemannian arclength t. 


2.6 Convolution and Haar-measure on SE{d ) 

In general a convolution of data U : SE(d) —> M with 
kernel K : SE{d) M is given by 

(K * U)(g) = f K(h~ 1 g) U(h ) d p(h) = 

SE(d ) 

/ / A((R')- 1 (x-x'),(R / )- 1 R)dx'd# so( d)(R / ), ( 25 ) 

R d SO(d) 

with d p{h) = dx'dpso(d) (R* , )> 

for all h = (x^R 7 ) E SE(d ), where Haar measure ~jJ 
is the direct product of the usual Lebesgue measure on 
with the Haar measure on SO(d). 


2.7 Gaussian Smoothing and Gradient on SE{d ) 

We define the regularized data 

V := G s * U, (26) 

where s = (s p ,s Q ) are the spatial and angular scales 
respectively of the separable Gaussian smoothing kernel 
defined by 

G s (x, R) := G?‘(x) Gf o d_1 (Ra). (27) 

This smoothing kernel is a product of the heat kernel 

ii X ii 2 

G^ d p (x) — s 4S y /2 on 1R. d centered at 0 with spatial scale 

s p > 0, and a heat kernel Gf d *(Ra) on centered 
around a E S d_1 with angular scale s Q > 0. 
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By definition the gradient VIZ of image data U : 
SE{d ) -A M is the Riesz representation vector of the 
derivative d U : 


W := % x dU = E fjL- 2 AiUAi-\ 

= M lt -,{A 1 U,...,A nd U) T , 


n d 

- E AjUAi 

j=d+1 


relying on as defined in (24). Here, following stan¬ 
dard conventions in differential geometry, 0” 1 denotes 
the inverse of the linear map associated to the metric 
tensor (23). Then, the Gaussian gradient is defined by 

W S U := VF = V(G S *U) = VG S * U. (29) 


2.8 Horizontal Exponential Curves in SE(d ) 



Fig. 4: Left: horizontal exponential curve 7 ^ in SE( 2) 
with c = (1, 0,1). Its projection on the ground plane re¬ 
flects co-circularity, and the curve can be obtained by 
a lift (30) from its spatial projection. Right: the distri¬ 
bution A of horizontal tangent vector fields as a sub¬ 
bundle in the tangent bundle T(SE( 2 )). 


Typically, in the distribution U (e.g. if U is an orien¬ 
tation score of a grey-scale image) the mass is concen¬ 
trated around so-called horizontal exponential curves 
in SE(d) (see Fig. [3|. Next we explain this notion of 
horizontal exponential curves. 

A curve t (x(t),y(t)) G M 2 can be lifted to a 
curve 1 1 —y 7 (t) = (x(t),y(t), 6 (t)) in SE( 2 ) via 

0(t) = arg {x(t) + i y{t)}. (30) 

Generalizing to d > 2, one can lift a curve t i-A x(t) G 
towards a curve t i-A 7 (t) = (x(£), n(£)) in M. d x S d ~ 1 
by setting 

n(t) = ||x(t)||- 1 x(t). 

A curve t i-A x(t) can be lifted towards a family of lifted 
curves t t-A 7 (t) = (x(t), R n (t)) into the roto-translation 
group SE(d) by setting R n (t) G SO(d ) such that it 
maps reference axis a onto n (t): 

Rn( t) a = n(t) = ||x(/) r'x(/). (31) 


— For d = 3, we impose the constraint: 


7 (t) G with A := span^ {A 3 , A 4 , A}, (33) 


where ^4 3 = n • Vrs, since then spatial transport is 


always along n which is required for for (31). 


Curves 7 (t) satisfying the constraint (32) for d = 2 , and 
( |33| ) for d = 3 are called horizontal curves. Note that 
dim(Z\) = d. 

Next we study how the restriction applies to the 
particular case of exponential curves on SE(d). 

— For d = 2 horizontal exponential curves are obtained 
from (18), (19), (20) by setting c 2 = 0. 

— For d = 3, we use a different reference axis a, and 


horizontal exponential curves are obtained from (16) 
by setting c 1 = c 2 = c 6 = 0 . 


If exponential curves are not horizontal, then we in¬ 
dicate how much the local tangent of the exponential 
curve points outside the spatial part of A, by a ‘de¬ 
viation from horizontality angle’ x? which is given by: 


Here we use R n to denote any rotation that maps refer¬ 
ence axis a onto n. Clearly, the choice of rotation is not 
unique for d > 2, e.g. if d = 3 then R n R a5Q! a = a re¬ 
gardless the value of <a, where R a ,« denotes the counter¬ 
clockwise 3D rotation about axis a by angle a. 


Next we study the implication of restriction (31) on 
the tangent bundle of SE(d). 

— For d = 2, we have restriction x(t) = (x(t),y(t)) = 
||x(t)||(cos 6 > (t), sind(t)), i.e. 


7 G A |~ , with A = spanjcos#c^ + sin 0d y ,de} 
= span{^i,^4 3 }, 


(32) 


where A denotes the so-called horizontal part of tan¬ 
gent bundle T(SE( 2 )). See Fig.[4j 


X = arccos 


-(i) 


|c«|| 


(34) 


Example 4 In case d = 2 we have 77,2 = 3, a = (1, 0) T . 
The horizontal part of the tangent bundle A is given 


by (32), and horizontal exponential curves are obtained 
from (18) by setting c 2 = 0. For exponential curves in 


general, we have deviation from horizontality angle 


X = arccos 


vW 


r.212 


(35) 


An exponential curve in SE( 2) is horizontal if and only 
if X = 0* See Fig. [4j where in the left we have depicted 
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a horizontal exponential curve and where in the right 
we have visualized distribution A. 


tion between the normalized gauge frame and the left- 
invariant vector field frame is given by 


Example 5 In case d = 3, we have ns = 6 , a = (0, 0,1) T . 
The horizontal part of the tangent bundle is given by 


(33), and horizontal exponential curves are character¬ 


ized by c 3 , c 4 , c 5 whereas c 1 = c 2 = c 6 = 0. By Eq. (22) 


these curves have zero torsion |r| = 0 and constant cur- 

/( c 4\2_^_f c 5\2 

vature ——^ 3 -^-and thus they are planar circles. For 

exponential curves in general, we have deviation from 
horizontality angle 


X = arccos 


>/fc 


112 


+ |c 2 | 2 + |c ; 


,312 


An exponential curve in SE( 3) is horizontal if and only 
if x = 0 and c 6 = 0 . 


3 From Exponential Curve Fits to Gauge 
Frames on SE(d ) 


B ■■= (R c ) t m; 1 a 


(36) 


with A := (Al,^I. 2 ,- 43 ) T , B := (£>i,S 2 ,S 3 ) t , and with 
rotation matrix 


R c = R 2 R 1 E SO( 3), with 

/cosx — sinx 0\ / cos v 0 sin v N 

R 2 = [ sin x cos x 0 ] , Ri = I 0 10 

\ 0 0 1/ \ — sin v 0 cos v 1 


(37) 


where the rotation angles are the deviation from hori¬ 
zontality angle x and the spherical angle 


v = arcsm 


e [-71-/2,71-/2]. 


Recall that x is given by (35). The multiplication M “ 1 A 
ensures that each of the vector fields in the locally adap¬ 
tive frame is normalized w.r.t. the 0 ^-metric, recall 

( 23 k. 


In Section [5] and Section [ 6 ] we will discuss techniques 
to find an exponential curve 7 g (t) that fits the data 
U : SE{d) R locally. Let c (g) = ( 7 ^)'( 0 ) be its 
tangent vector at g. 

In this section we assume that the tangent vector 
c(g) = (c^^(^), c^ 2 \g)) T E R d+rd = M nd is given. 
From this vector we will construct a locally adaptive 
frame {B\\ g ,..., B nd \ g }, orthonormal w.r.t. 0 ^-metric 
in such a way that: 

1 . the main spatial generator (Ai for d = 2 and Ad for 

d> 2 ) is mapped onto Bi\ g = ° l {g) Ai\ g , 

i= 1 

2 . the spatial generators {Bi | }f =2 are obtained from 

the other left-invariant spatial generators {Ai\ g }f =1 
by a planar rotation of a onto by angle x* I n 

particular, if x = 0 , the other spatial generators do 
not change their direction. This allows us to still dis¬ 
tinguish spatial generators and angular generators 
in our adapted frame. 

Next we provide for each g E SE(d) the explicit con¬ 
struction of a rotation matrix R c( ^ and a scaling by 
M^-i on T g (SE(d)), which maps frame {A\\ g ,..., A nd \ g } 
onto {Bi\ g ,...,B nd \ g }. 

The construction for d > 2 is technical and provided 
in Theorem [A| in Appendix [A[ However, the whole con¬ 
struction of the rotation matrix R c via a concatenation 
of two subsequent rotations is similar to the case d = 2 
that we will explain next. 

Consider d — 2 where the frames {Ai, A 2 , As} and 
{B \, jB 2 , Bs} are depicted in Fig. [ 5 ] The explicit rela- 


Remark 2 When imposing isotropy (w.r.t. the metric 
0 M ) in the plane orthogonal to £> 1 , there is a unique 
choice R c mapping (1,0, 0) T onto (/zc 1 , /m 2 , c 3 ) T such 
that it keeps the other spatial generator in the spa¬ 
tial subspace of T g (SE( 2)) (and with x = 0 S 2 = 
/i _ 1 7l 2 ). This choice is given by (37). 



Fig. 5: Locally adaptive frame {Bi\ , S 2 |^ , Bs\ g } (in 
blue) in T g (SE( 2)) (with g placed at the origin) is 
obtained from frame {Ai\ g , A 2 \ g , As\ g } (in red) and 
c(g), via normalization and two subsequent rotations 
R c = R 2 Ri, see Eq. ( |36| ) , revealing deviation from hor¬ 
izontality x in i?i, spherical angle v in Eq. (37). Vec¬ 
tor field A\ takes a spatial derivative in direction n, 
whereas B\ takes a derivative along the tangent c of 
the local exponential curve fit. 
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The generalization to the d-dimensional case of the con¬ 
struction of a locally adaptive frame from {AiYi=\ 

and the tangent vector c of a given exponential curve fit 
7 g(-) to data U : SE(d) -A M is explained in Theorem [7] 
in Appendix [Aj 


4 Exponential Curve Fits in 

In this section we reformulate the classical construc¬ 
tion of a locally adaptive frame to image / at location 
x G M d , in a group-theoretical way. This reformulation 
seems technical at first sight, but helps in understand¬ 
ing the formulation of projected exponential curve fits 
in the higher dimensional Lie group SE(d). 


4.1 Exponential Curve Fits in of the 1 st Order 


We will take the structure tensor approach [THUHT)] , which 
will be shown to yield first-order exponential curve fits. 
The Gaussian gradient 

V s / = VG S * /, (38) 

with Gaussian kernel 

G s (x) = , (39) 

is used in the definition of the structure matrix: 


s s ’ p (f) = G p * V s /(V s /) 3 


(40) 


with s = \o 2 , and p = \a 2 p the scale of regularization 
typically yielding a non-degenerate and positive definite 
matrix. In the remainder we use short notation S s,p := 
S s,p (/). The structure matrix appears in solving the 
following optimization problem where for all x G l d we 
aim to find optimal tangent vector 


c*(x) = argmin f G p (x — x')|V s /(x') • c 

c e R d , R d 
l|c||~i 

= argmin c T S s,p (x)c. 

c G 

Hell = 1 


2 dx' 


(41) 


In this optimization problem we find the tangent c* (x) 
which minimizes a (Gaussian) weighted average of the 
squared directional derivative | V s /(x')-c | 2 in the neigh¬ 
borhood of x. The second identity in (41), which di¬ 


rectly follows from the definition of the structure ma¬ 


trix, allows us to solve optimization problem (41) via 
the Euler-Lagrange equation 


S s ’ p (x ) c*(x) = Aic*(x), (42) 

since the minimizer is found as the eigenvector c*(x) 
with the smallest eigenvalue Ai. 


Now let us put Eq. (41) in group-theoretical form by 
reformulating it as an exponential curve fitting prob¬ 
lem. This is helpful in our subsequent generalizations 
to SE(d). On exponential curves are straight lines: 


7 xW = x + exp Rd (tc) = x + tc, 


(43) 


and on T(M d ) we impose the standard flat metric tensor 
0(c, d) = Yli=i ^d l . In (41) we replace the directional 
derivative by a time derivative (at t = 0 ) when moving 
over an exponential curve: 


c*(x) = argmin 


ceM d ,|| 

c\\ = l 

9 

/ G p (x - x') 

|(G 5 */)( 7 ^,xW )| f=0 

dx', 

R d 



(44) 


where 


t •-> 7x',xW = 7xW - x + x ' = 7x'(h (45) 

Because in ( |41| we average over directional derivatives 
in the neighborhood of x we now average the time 
derivatives over a family of neighboring exponential curves 
7 x 7 xM? which are defined to start at neighboring posi¬ 
tions x' but having the same spatial velocity as 7 £(£). 

In R d the distinction between 7 ^, x (t) and 7 x ,(t) is not 
important but it will be in the SE(d)- case. 

Definition 1 Let c*(x) G T x (M d ) be the minimizer in 
( |44| ). We say 7 x (t) = x + exp Rd (tc*(x)) is the first- 
order exponential curve fit to image data / : -A M 

at location x. 


4.2 Exponential Curve Fits in R d of the 2 nd Order 


For second-order exponential curve fits we need the 
Hessian matrix defined by 

(H s (/))( x ) = [d Xj d Xi (G s * /)(x)] , (46) 

with G s the Gaussian kernel given in Eq. (|39|). From now 
on we use short notation H s := H s (/). When using the 
Hessian matrix for curve fitting we aim to solve 

c*(x) = argmin |c T H s (x)c|. 

c e (47) 

M = 1 


In this optimization problem we find the tangent c* (x) 
which minimizes the second-order directional derivative 
of (Gaussian) regularized data G s *f. When all Hessian 
eigenvalues have the same sign we can solve the opti¬ 
mization problem (47) via the Euler-Lagrange equation 


H s ( x ) c*(x) = Aic*(x), 


(48) 
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and the minimizer is found as the eigenvector c*(x) 
with the smallest eigenvalue Ai. 

Now, we can again put Eq. (47) in group-theoretical 
form by reformulating it as an exponential curve fitting 
problem. This is helpful in our subsequent generaliza¬ 
tions to SE(d). We again rely on exponential curves 
as defined in (|43|). In (|47|) we replace the second order 


directional derivative by a second order time derivative 
(at t = 0 ) when moving over an exponential curve: 


ones used in Section [4] we will, for all three cases, first 
pose an optimization problem, and then give its solu¬ 
tion in a subsequent theorem. 

In this section we rely on the group-theoretical tools 
explained in Section [ 2 ] (only the case d= 2 ), listed in 
subtables E.l and E.2 in our table of notations. Fur¬ 
thermore we introduce notations listed in the first part 
of subtable E.3. 


c*(x) = argmin 

ceM d ,||c||=i 


^2 (G s * /)(7xW) 


t =o 


(49) 


Remark 3 In general the eigenvalues of Hessian matrix 
H s do not have the same sign. In this case we still take 
c*(g) as the eigenvector with smallest absolute eigen¬ 
value (representing minimal absolute principal curva¬ 
ture), though this no longer solves (47). 


Definition 2 Let c*(x) G T x (M d ) be the minimizer in 
( |49| ). We say y x (t) = x + exp R d(tc*(x)) is the second- 
order exponential curve fit to image data / : R d M 
at location x. 


5.1 Neighboring Exponential Curves in SE(2 ) 


Akin to (45) we fix reference point g G SE( 2) and ve¬ 
locity components c = c(g) G M 3 , and we shall rely on a 
family { 7 ^ g } of neighboring exponential curves around 
7 g. As we will show in subsequent Lemma [l] neighbor¬ 
ing curve 7 ^ departs from h and has the same spatial 
and rotational velocity as the curve 7 ^ departing from 
g. This geometric idea is visualized in Fig.[ 6 j where it is 
intuitively explained why one needs the initial velocity 
vector R/j-ipC, instead of c in the following definition 
for the exponential curve departing from a neighboring 
point h close to g. 


Remark 4 In order to connect optimization problem 


(49) with the first order optimization (44) we note that 


(49) can also be written as an optimization over a fam¬ 


ily of curves 7 X , x defined in (45): 


*(x) = 


arg min 

c G 

ll c ll = 1 


jG s (x-x')|(/)fc jX (t)) 

R d 


t =0 


dx'. 


Definition 3 Let g G SE(2) and c = c (g) G M 3 be 
given. Then we define the family { 7 ^ g } of neighboring 
exponential curves 


* %g(t) ■= lh 


R, _i c . . 
h 9 (*), 


( 51 ) 


with rotation-matrix R^-i^ G SO(3) defined by 


1 g 


^(R') r R 0^ 


(52) 


because of linearity of the second-order time derivative. 

5 Exponential Curve Fits in SE( 2) 


for all g = (x, R) G SE{2) and all h = (x',R 7 ) G 
SE( 2), with R, R 7 G SO(2) a counterclockwise rotation 
by respectively angle 0 and 0 '. 


As mentioned in the introduction we distinguish be¬ 
tween two approaches: a first order optimization ap¬ 
proach based on a structure tensor on SE( 2 ), and a 
second order optimization approach based on the Hes¬ 
sian on SE( 2 ). The first order approach is new while 
the second order approach formalizes the results in [351 
129] . They also serve as an introduction to the new, more 
technical, 5'E(3)-extensions in Section [ 6 j 

All curve optimization problems are based on the 
idea that a curve (or a family of curves) fits the data 
well if a certain quantity is preserved along the curve. 
This preserved quantity is the data U (q(t)) for the first 
order optimization, and the time derivative ^ 1 /( 7 (t)) 
or the gradient V£7 (7 (t)) for the second order optimiza¬ 
tion. After introducing a family of curves similar to the 


Lemma 1 Exponential curve 7 ^ departing from h G 
SE( 2 ) given by (51) has the same spatial and angular 
velocity as exponential curve 7 ^ departing from g G SE{2). 

On the Lie algebra level; we have that the initial 
velocity component vectors of the curves 7 ^ and 7 ^ 
relate via c R h -i g c. 

On the Lie group level; we have that the curves them- 


selves jg(-) = (x g (-),R g (-)), 7 h,gO) = ( x h(-),Rh(-)) re- 
late via 


•Bh(f) - X g(f) X T ^ 1 /ro\ 

R h (t) = RglfyR-'R! 9 h (t) = 9 g {t ) -9 + 6'. [06> 

Proof The proof follows from the proof of a more gen¬ 
eral theorem on the SE( 3) case which follows later (in 
Lemma [3| . 
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Remark 5 Eq. (53) is the extension of Eq. (45) on M 2 to 


the SE( 2 ) group. 

Additional geometric background is given in Appendix[B] 




Fig. 6 : Family of neighboring exponential curves, given 
a fixed point g G SE( 2) and a fixed tangent vector 
c = c (g) G T g (SE( 2)). Left: Our choice of family of ex¬ 
ponential curves 7 ^ g for neighboring h G SE( 2). Right: 
Exponential curves 7 ^ with c = c(g) are not suited for 
local averaging in our curve fits. The red curves start 
from g (indicated with a dot), the blue curves from 
h^g. 


5.2 Exponential Curve Fits in SE( 2 ) of the 1 st Order 
For first-order exponential curve fits we solve an opti¬ 


mization problem similar to (44) given by 


( 54 ) 

with V = G s * C/, g = (x, R), h = (x', R/) and dp(h) = 
dx / d/i 5 o( 2 )(R' / ) = dx'd O'. Here we first regularize the 
data with spatial and angular scale s = (s p ,s Q ) and 
then average over a family of curves where we use spa¬ 
tial and angular scale p = (p p ,p G ). Here s p , p p > 0 are 
isotropic scales on M 2 and s 0 ,p 0 > 0 are scales on S 1 
of separable Gaussian kernels, recall ([27]). Recall also 
(24) for the definition of the norm || • || M . When solv¬ 


ing this optimization problem the following structure 
matrix appears 


(S s ’ p (U))(g) = f G p {h-'g). 

SE{ 2) 


~ T 


(55) 


R h - lq VV(h)(yV(h)) T R h -i g dp(h). 


In the remainder we use short notation S s,p := S S,P (U). 
We assume that C7, p, s, and p, are chosen such that 
S s,p (g) is a non-degenerate matrix. The optimization 
problem is solved in the next theorem. 

Theorem 1 (First Order Fit via Structure Tensor) 

The normalized eigenvector M p c*(g) with smallest eigen¬ 
value of the rescaled structure matrix 
M ll S s,p (g)M /Jj provides the solution c*(g) to optimiza¬ 
tion problem ( 54 ] ). 

Proof We will apply four steps. In the first step we 
write the time-derivative as a directional derivative, in 
the second step we express the directional derivative in 
the gradient. In the third step we put the integrand in 
matrix-vector form. In the final step we express our op¬ 
timization functional in the structure tensor and solve 
the Euler-Lagrange equations. 


For the first step we use (51) and the fundamen- 


application of (13): 


tal property (17) of exponential curves such that via 


h,g 


t =0 


( 0 ),4g i9 (o)> 


(dy^R^-igC) 


(56) 


where we use short notation R^-i^c = 5^ i _ 1 (R/ l -i p c)i4i|/ l . 
In the second step we use the definition of the gra¬ 


dient (28) and the metric tensor (23) to rewrite this 


expression to 

(dV'| /l ,R h -i g c)) = (VV(h),R h -i g c) 


(57) 


Then, in the third step we write this in vector- 
matrix form and obtain 


®/*lh(W(A), Rft-ipc) 


c r M^R l lg VV(h) 


c ( g) = argmin [G p (h 1 g) 

cGM 3 , J 


t =0 

2 

dp{h), 

I|c|| m = i se ^ 





= c T M^Rl_ lg VV(h)(VV(h)) T R h -i g M, 2 C, 


(58) 


where we used the fact that M ^2 and R^_i^ commute. 


Finally, we use the structure tensor definition (55) 


to rewrite the convex optimization functional in (54) as 


t =0 


dp(h) 


(59) 


f(c) := f G p {h~ 1 g) &V(% g (t)) 

SE( 2) 

= c T M |[t 2S s ’ p M /i 2 C , 

while the boundary condition Hc^ = 1 can be written 
as 


<p(c) := c t M m 2 C — 1 = 0. 


(60) 


The Euler-Lagrange equation reads V£(c*) = AiV<p(c*), 
with Ai the smallest eigenvalue of M Ai S s,p (^ f )M /i and 
we have 


M^S^M^c*^) | Ax M M 2 C *(g) 
M^ p (g)M p (M^c*(g)) = A^M^c*^)) 


(61) 
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from which the result follows. 

The next remark explains the frequent presence of the 


□ Remark 8 Optimization problem (63) can also be writ- 


M w matrices in (69) 


Remark 6 The diagonal matrices enter the func¬ 
tional due to the gradient definition (28), and they en¬ 
ter the boundary condition via ||c||^ = c^M^c = 1. In 


both cases they come from the metric tensor (23). Pa¬ 
rameter fi which controls the stiffness of the exponential 
curves has physical dimension [Length] -1 . As a result, 
the normalized eigenvector M m c *(g) is, in contrast to 
c*(g), dimensionless. 


5.3 Exponential Curve Fits in SE( 2) of the 2nd Order 

We now discuss the second order optimization approach 
based on the Hessian matrix. At each g E SE( 2) we 
define a 3 x 3 non-symmetric Hessian matrix 

(H s ([7))(2) = [AjMVKg)] , with V = G a *U, (62) 

and where i denotes the row index and where j denotes 
the column index, and with G s a Gaussian kernel with 


isotropic spatial part as described in Eq. (27). In the 
remainder we write H s := H S (U). 

Remark 1 As the left-invariant vector fields are non- 
commutative there are many ways to define the Hes¬ 
sian matrix on SE( 2), since the ordering of the left- 
invariant derivatives matters. From a differential geo¬ 


metrical point of view our choice (62) is correct, as we 
motivate in Appendix |D| 

For second-order exponential curve fits we consider 2 
different optimization problems. In the first case we 
minimize the second order derivative along the expo¬ 
nential curve: 



d? ~ 



c*(g) = arg min 

cGM 3 ,||c|| Al = l 

^ v m)) 

t =0 



(63) 


In the second case we minimize the norm of the first or¬ 
der derivative of the gradient of the neighboring family 
of exponential curves: 


c *(g) = arg min 

f G p (h *,/)• 


c€R3,||c|| m =1 

SE( 2) 

J Mh) ■ 

®m(IvC(7 Zjt)) 

t=0 >|VE(7^(f)) 


( 64 ) 


with again V = G s * U. 


ten as an optimization problem over the neighboring 
family of curves, as it is equivalent to problem: 


c*(g) = arg min 

CP 3 


c e 




d‘ 


GsihT^g) —U{% g {t)) 


= 1 SE{ 2) 


dp(h) 


t =o 


(65) 


In the next two theorems we solve these optimiza¬ 
tion problems. 

Theorem 2 (Second Order Fit via Symmetric Sum 
Hessian) Let g E SE( 2) be such that the eigenvalues 
of the rescaled symmetrized Hessian 

l -M;\H s {g) + {H s {g)) T )M- 1 

have the same sign. Then the normalized eigenvector 
M fJj c*(g) with smallest eigenvalue of the rescaled sym¬ 
metrized Hessian matrix provides the solution c*(g ) of 


optimization problem (63). 


Proof Similar to the proof of Theorem [l] we first write 
the time derivative as a directional derivative using Eq. 


(13). Since now we have a second order derivative this 


step is applied twice: 

&nr g (t)) 


t =o 


Tt £ C *A^(7 g(t)) 
2=1 


t =0 


E Aj{AV){g) 

i,j =1 


( 66 ) 


Then we write the result in matrix-vector form and split 
the matrix in a symmetric and anti-symmetric part 


E Aj{AV){g) 

i,j =1 


|c t H s ( 5 )c| 


= |ic r (H s ( 5 ) + (H s ( 5 ))V 

+ |c t (H s ( 5 ) - (U s (g)) T ) 
= h |c t (H s (5() + (H s ( 5 )) t )c| , 


(67) 


where only the symmetric part remains. Finally, the 
optimization functional in ( [63] ) (which is convex if the 
eigenvalues have the same sign) can be written as 


£(c) : = famm 


t =o 


Uc T (H s (g) + (H s (g)) T )c\. 


( 68 ) 


Again we have the boundary condition c) = c J M (l 2 C- 
1 = 0. The result follows using the Euler-Lagrange for¬ 
malism V£(c*) = AiV<y9(c*): 

|(H S ( 5 ) + (U s (g)) T )c*(g) = A, M^c*(g) * 

±M-pK s (g) + (H s (g)) T )M-i(M,c*(g)) 

= X 1 (M,c*(g)), 


( 69 ) 
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which boils down to finding the eigenvector with mini¬ 
mal absolute eigenvalue |Ai| which gives our result. □ 


Theorem 3 (Second Order Fit via Symmetric 
Product Hessian) Let p p: p Q , s p , s Q > 0. The normal¬ 
ized eigenvector M p c*(g) with smallest eigenvalue of 
matrix 


M~P f G p (h~ 1 g) ■ R^-i g (H s (h)) T 

SE( 2) 

MpH°(h)R h -, g dJL(h) M~P 


(70) 


provides the solution c*(g) of optimization problem [6J±). 


Proof First we use the definition of the gradient (28) 
and then we again rewrite the time-derivative as a di¬ 
rectional derivative: 


ivv(r h pt)) _ = iE Av( 7 ^)K 2 Ab 


dt ^ JEi ^It git) 

t —u i~ i 

= c 3 AjAiV(h)npAi\ h 
i,3 =1 


t=0 


(71) 


for c = Rft-i^c, recall (52), and where pi = p for i = 
1, 2 and pi = 1 for i = 3. Here we use 7/^(0) = ft, and 
the formula for left-invariant vector fields ([l0| ). Now 
insertion of into the metric tensor ( |23| ) yields 

= c t (H s (/i)) t M“^H s (/i)c 
= c T R.7 l 5 (H s (/i)) T M“ 2 H s (/i)R^-i s c. 


(72) 


Finally, the convex optimization functional in (64) can 
be written as 


£(c) := c T 


f G p (h- 1 g)-R h - lg (H'(h)) T 

SE(2) (73) 


MpH s (h)R h -i g djl(h) c. 


Again we have the boundary condition y>(c) = c T M /t 2 C— 
1 = 0 and the result follows by application of the Euler- 
Lagrange formalism: V£(c*) = AiV(/?(c*). □ 


6 Exponential Curve Fits in SE( 3) 

In this section we generalize the exponential curve fit 
theory from the preceding chapter on SE( 2) to SE(3). 
Because our data on the group SE( 3) was obtained 
from data on the quotient M 3 xi S 2 we will also discuss 
projections of exponential curve fits on the quotient. 

We start in Subsection |6.1| with some prelimenaries 
on the quotient structure (ph). Here we will also intro¬ 
duce the concept of projected exponential curve fits. 


Subsequently, in Subsection |6.2[ we provide basic the¬ 
ory on how to obtain the appropriate family of neigh¬ 
boring exponential curves. More details can be found in 
Appendix [B] In Subsection |6.3| we formulate exponen¬ 
tial curve fits of the first order as a variational prob¬ 
lem. For that we define the structure tensor on SE(3), 
which we use to solve the variational problem in Theo¬ 
rems [4] and |5j Then we present the two-fold algorithm 
for achieving torsion-free exponential curve fits. In Sub¬ 
section |6.4| we formulate exponential curve fits of the 
second order as a variational problem. Then we define 
the Hessian tensor on SE{3 ), which we use to solve the 
variational problem in Theorem [6j Again torsion-free 
exponential curve fits are accomplished via a two-fold 
algorithm. 

Throughout this section we will rely on the differen¬ 
tial geometrical tools of Section[2j listed in subtables [E|l 
and @2 in Appendix |EJ We also generalize concepts on 
exponential curve fits introduced in the previous sec¬ 
tion to the case d = 3 (requiring additional notation). 
They are listed in subtable |E|3 in Appendix [EJ 


6.1 Preliminaries on the quotient M 3 xi S 2 . 

Now let us set d = 3, and let us assume input U is given 
and let us first concentrate on its domain. This domain 
equals the joint space M 3 xi S 2 of positions and ori¬ 
entations of dimension 5, which we identified with a 5- 
dimensional group quotient of SE( 3), where SE( 3) is of 
dimension 6 (recall (|3|). For including a notion of align¬ 
ment it is crucial to include the non-commutative rela¬ 
tion in between rotations and translation, and not to 
consider the space of positions and orientations as a hat 
Cartesian product. Therefore we model the joint space 
of positions and orientations as the Lie group quotient 
(J3|, where 

50( 2) = Stab (a) = {R G 50(3) | Ra = a} 

for reference axis a = e z = (0, 0,1) T . Within this quo¬ 
tient structure two rigid body motions g = (x, R),g' = 
(x',R') G SE( 3) are equivalent if 

g' ~ g W)~ l g e { 0 } x 50(2) 

x-x'= 0 and aeSO(2 ) : (R') _1 R = R e *,a- 

Furthermore, one has the action 0 of g — (x, R) G 
SE( 3) onto (y, n ) G M 3 x S 2 , which is defined by 

9 © (y, n ) = ( x , R) © (y, n) := (x + Ry, Rn). (74) 

As a result we have 

g' ~ 9 O g' © (0, a) = g © (0, a). 
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Thereby, a single element in M 3 x S 2 can be consid¬ 
ered as the equivalence class of all rigid body motions 
that map reference position and orientation (0, a) onto 
(x, n). Similar to the common identification of S 2 = 
SO(3)/SO(2), we denote elements of the Lie group quo¬ 
tient M 3 x S 2 by (x, n). 


6 .1.1 Legal Operators 

Let us recall from Section [3] that exponential curve 
fits induce gauge frames. Note that both the induced 
gauge frame {Si,..., Bq] and the non-adaptive frame 
{Al, . .. ,Aq} are defined on the Lie group SE( 3), and 
cannot be defined on the quotient. Nevertheless, com¬ 
binations of them can be well-defined on M 3 x S 2 (e.g. 
A r s = A 2 + A 2 + A 2 is well-defined on the quotient). 
This brings us to the definition of so-called legal oper¬ 
ators, as shown in |30j Thm.l]. In short, the operator 
U i —y @(U) is legal (left-invariant and well-defined on 
the quotient) if and only if 

<p = <p o 7 Z ha for all a G [0, 2tt). 
o C g = C g o 3 for all g G SE( 3), 

recall (|7|), where 

ha * = (0,R Gz?a ). 

with the R e ^,a the counterclockwise rotation about e z . 
Such legal operators relate one-to-one to operators <L> : 
L 2 (M 3 x S 2 ) L 2 (M 3 x S 2 ) via 

U^$(U) ^ &*->$(&) =${U), 
relying consequently on §. 


(75) 

(76) 


6.1.2 Projected Exponential Curve Fits 


Action (74) allows us to map a curve y(-) = (x(-),R(-)) 
in SE( 3) onto a curve (x(-),n(-)) on M 3 x S 2 via 

(x(t), n(t)) := 7 (f) © (0, e z ) = (x(t), R (t) e z ). (77) 

This can be done with exponential curve fits 7 g~ c { ' J) (t) 
to define projected exponential curve fits. 

Definition 4 We define for g = (x, R n ) the projected 
exponential curve fit 

7*x,„)(i):=7f (s) (i)®(0,e,). (78) 

Lemma 2 The projected exponential curve jit is well- 
defined on the quotient, i.e. the right-hand side of (78) 
is independent of the choice of R n s.t. R n e z = n, if the 
optimal tangent found in our fitting procedure satisfies: 


(79) 


and for all g G SE( 3), with 

Z a ■■= R ° e J € SO( 6 ). (80) 

Proof For well-posed projected exponential curve fits 
we need the right-hand side of (78) to be independent 
of R n s.t. R n e z = n i.e. it should be invariant under 
g —> gh a • Therefore we have the following constraint on 
the fitted curves: 


7 f (9) (t) © (0,e z ) = 7 f h [ 9ha \t) © (0,e z ). 


(81) 


Then the constraint on the optimal tangent (79) follows 
from fundamental identity 


i.l C qha {-))=la aC {-)h, 


(82) 


which hold^] for all h a . We apply this identity ( [82] ) to 
the right-hand side of ( |8Tj ) and use the definition of O 
defined in ([74]) yielding: 


7 g i9 \t) © ( 0 , e z ) = 7 g“ c {9ha \t)h a © ( 0 , e z ) 

7 g i9 \t) © (0,e 2 ) = 7 g“ c (gh -\t) © (0, e 2 ) (83) 

t 

c *(g) = Z a c *{gh a ). 


Finally our constraint (|79j) follows from = Zq, 1 . □ 


6.2 Neighboring Exponential Curves in SE( 3) 


Here we generalize the concept of family of neighboring 
exponential curves (45) in the M d -case, and Definition [ 3 ] 
in the SE( 2)-case, to the SE( 3)-case. 


Definition 5 Given a fixed reference point g G SE( 3) 
and velocity component c = c (g) = (c^ (g), (g)) G 

M 6 , we define the family { 7 ^ ^(-)| of neighboring expo¬ 
nential curves by 




(84) 


with rotation matrix R^-i g G SO (6) defined by 


R^/i — 1 g 


/(R') t R 0 \ 

V 0 (r') t r7 


(85) 


for all g = (x, R), h = (x', R') € SE( 3). 


The next lemma motivates our specific choice of neigh¬ 
boring exponential curves. The geometric idea is visu¬ 
alized in Fig. [7] and is in accordance with Fig. [6] on the 
SE{2) case. 


Eq. (82) follows from (122) in App. B 


by setting Q = Z a . 


c*(gh a ) = Zlc*{g), for all a e [0,2?r], 
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Fig. 7: Illustration of the family of curves 7 ^ in SE( 3). 
Left: The (spatially projected) exp-curve t P^jgit), 
with g = (x, R n ) in red. The frames indicate the rota¬ 
tion part Pso(3)lg{t), which for clarity we depicted only 
at two time instances t. Middle: neighboring exp-curve 
t i-A 7 g^(t) with h = (x', R n ), x ^ x' in blue, i.e. neigh¬ 
boring exp-curve departing with same orientation and 
different position. Right: exp-curve t 7 g h (t) with 

h = (x, R n /), n' 7 ^ n, i.e. the neighboring exp-curve 
departing with same position and different orientation. 


Lemma 3 Exponential curv e 7 ? g departing from h = 
(od , R') G SE{3) given by ( 84) has the same spatial 
and rotational velocity as exponential curve 7 ^ depart¬ 
ing from g = (x, R ) G SE( 3). 

On the Lie algebra level; we have that the initial 
velocity component vectors of the curves 7 ^ and 7 ^ 
relate via c i-» R h -i g c. 

On the Lie group level; we have that the curves them¬ 


selves Jg(-) = (x g {-),Rg(-)), 7 S (-) = (x h (-),R h {-)) re- 
late via 


x hif) — x gif) x-\-x /, . ^ 

R h (t) = R g (t)R~ % R'. 1 j 

Proof See Appendix [Bj 

Remark 9 Lemma [ 3 ] extends Lemma [l] to the SE( 3) 
case. When projecting the curves 7 ^ and 7 ^ into the 
quotient, one has that curves 7 ^ 0 ( 0 , a), and 7 ^ 0 ( 0 , a) 
in M 3 xi S 2 carry the same spatial and angular velocity. 


6.3 Exponential Curve Fits in SE( 3) of the 1st Order 

Now let us generalize the first-order exponential curve 
fits of Theorem[l]to the setting of M 3 xi S' 2 . Here we first 
consider the following optimization problem on SE( 3) 
(generalizing (|44|)): 


c*(g) = arg min G p (h x g) 

ceM 6 , J 


2 

t =0 

d p(h), 

||o|U = 1,^(3) 




0 

01 

II 

0 





( 88 ) 


Recall that || • || M was defined in (24), V in (26) and fi 
in (25). The reason for including the condition c 6 = 0 
will become clear after defining the structure matrix. 


6.3.1 The Structure Tensor on SE( 3) 


We define structure matrices S s,p of U by 
(s s ’ p (U))(g) = f Gpih-ty. 

SE( 3) 

RE Q V s U(h)(V s U(h)) T R k - lg dJl(h), 


(89) 


where we use matrix R^-i g defined in Eq. (85). Again 
we use short notation S S)P := S S,P (U). 

Remark 11 By construction 0 and ( flQ| ) we have 


so the null space of our structure-matrix includes 

AT := span{(0, 0, 0, 0, 0,1) T }. (90) 

Remark 12 We assume that s = (s p , s Q ) and function U 
are chosen in such a way that the null space of the struc¬ 
ture matrix is precisely equal to A f (and not larger). 

Due to the assumption in Remark [12] we need to impose 
the condition 


Remark 10 In order to construct the family of neigh¬ 
boring exponential curves in SE( 3) one applies the trans¬ 
formation c Rft-i^c in the Lie algebra. Such a trans¬ 
formation preserves the left-invariant metric: 


1 = 0 


3oP) 


(% c w,% c w) = ®7 (t) (7 c h jt),r h jt)), 

(87) 


for all h G SE( 3) and all t G M. For further differential 
geometrical details see Appendix [B] 


c 6 = o ^ 4o(o) n W = 0 


(91) 


in our exponential curve optimization to avoid non¬ 
uniqueness of solutions. To clarify this, we note that 


the optimization functional in (88) can be rewritten as 


£(c) := c T M^S s ’P(g)M^c, 


as we will show in the next theorem where we solve the 
optimization problem for first-order exponential curve 
fits. Indeed, for uniqueness we need ( |9lj ) as otherwise 
we would have £ (c + M m - 2 Co) = £(c) for all Cq G Af. 
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Theorem 4 (First Order Fit via Structure Ten¬ 
sor) The normalized eigenvector M /Jj c*(g) with small¬ 
est non-zero eigenvalue of the rescaled structure matrix 
M /d S s ’ p (g)M ld provides the solution c*(g) to optimiza¬ 
tion problem [88 


Proof All steps (except for the final step of this proof, 
where the additional constraint c 6 = 0 enters the prob¬ 
lem) are analogous to the proof of the first order method 
in the SE(2) case: the proof of Theorem [l] We will now 
shortly repeat these first steps. First we rewrite the 
time derivative as a directional derivative which is then 
rewritten to the gradient 


(V(7kjt))) 


t =0 


< d ^k 9 (0),7U°)> 

(dF|/i,R*-i ff c) 
®A(VV-(>0, R fc -i ff c) 


(92) 


We then put this result in matrix-vector form: 

= c T M^Rl_ lg (WV(h))(VV(h)) T R h -i g M^c. 

This again yields the following optimization functional 


(93) 


W =JsE { 3)G P (h- 1 g)\£V(r h , g (t)) 

= c T M^S s ’P(g)M^c. 


t=o 


dp(h) 


(94) 


So, just as in the SE( 2)-case we have the following 
Euler-Lagrange equations: 


M^S s ’ p (g)M^c*(g) = AxM^c *(g) 
M At S s ’ p ( 5 )M A1 (M At c*( 5 )) = A! (M M c*(<?)). 


(95) 


Again the second equality in ([95]) follows from the first 
by multiplication by M” 1 . 

Finally, the constraint c 6 = 0 is included in our op¬ 


timization problem (88) to excluded the null space (90) 


from the optimization, therefore we take the eigenvec¬ 
tor with the smallest non-zero eigenvalue providing us 
the final result. □ 

6.3.2 Projected Exponential Curve Fits in M 3 x S 2 

In reducing the problem to M 3 x S 2 we first note that 
S s ’P(gh a ) = ZlS s ’P(g)Z a , (96) 


with Z a defined in Eq. (80), and where we recall h a = 
(0,Re„a). 

In the following theorem we summarize the well- 
posedness of our projected curve fits on data V : M 3 x 
S 2 M and use the quotient structure to simplify the 
structure tensor. 


Theorem 5 (First Order Fit and Quotient Struc¬ 
ture) Let g = ( x , R n ) and h = (a/, R n >) where R n and 
R n > denote any rotation which maps e z onto n and n' 
respectively. Then, the structure tensor defined by 
can be expressed as 

s s ’P(g) = 27 rf f Gfpx- o') Gf 0 {R T n ,n) 

r 3 s 2 (97) 

R Z- lg VV(h) (\/V(h)) T R h -i g da(n')d^. 

The normalized eigenvector M^c*(x, R n ) with small¬ 
est non-zero eigenvalue of the rescaled structure matrix 


M^S s ' p {g)M^provides the solution of (88) and defines 
a projected curve fit in M 3 x S 2 ; 

) © (0, e z ), (98) 

which is independent of the choice of R n > and R n . 

Proof The proof consists of two parts. First we prove 


that (97) follows from the structure tensor defined in 


(89). Then we use Lemma[2]to prove that our projected 
exponential curve fit ( [98] ) is well-defined. For both we 
use Theorem [4] as our venture point. 

For the first part of the proof we note that the in¬ 


tegrand in the structure tensor definition Eq. (89) is 


invariant under h hh a = h( 0,R ez?a ) on the integra¬ 
tion variable. To show this we first note that Z a defined 
in (80), satisfies Z a (Z a ) T = I. Furthermore, we have 


S7V(hh a ) = ZaVV(h), R(fch a )-i g — 

and G p (hh a ) = G p (h). Therefore integration over third 
Euler-angle a is no longer needed in the definition of the 


structure tensor (89) as it just produces a constant 2i r 
factor. 

For the second part we apply Lemma [2] and thereby 
it remains to be shown that condition c*(gh a ) = Z^c *(g) 
is satisfied. This directly follows from ([Ml): 


(99) 


S 3lp {gh ot )c*(gh ot ) = X 1 c*(gh a ) 

t 

ZlS s ’P(g)Z a c*(gh a ) = A lC *(gh a ) 

S s ’ p (g) (Z a c*(gh a )) = Ai (Z Q c *(gh a )) 

t 

Z a c* (gh a ) = c*(g), 
which shows our condition. 


□ 


6.3.3 Torsion-free Exponential Curve Fits of the 1st 
Order via a Two-fold Approach 

Theorem [4] provides us exponential curve fits that pos¬ 
sibly carry torsion. From Eq. (22) we deduce that the 
torsion norm of such an exponential curve fit is given by 
l r l = || c (i)|| (c 1 c 4 +c 2 c 5 +c 3 c 6 )|^|. Together with the fact 
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that we exclude the null space A f from our optimiza¬ 
tion domain by including constraint c 6 = 0 , this results 
in insisting on zero torsion along horizontal exponen¬ 
tial curves where c 1 = c 2 =0. Along other exponential 
curves torsion appears if c x c 4 + c 2 c 5 ^ 0 . 

Now the problem is that insisting, a priori, on zero 
torsion for horizontal curves while allowing non-zero 
torsion for other curves is undesirable. On top of this, 
torsion is a higher order less-stable feature than curva¬ 
ture. Therefore we would like to exclude it altogether 
from our exponential curve fits presented in Theorem [4] 
and Theorem [5| by a different theory and algorithm. 
The results of the algorithm show that even if struc¬ 
tures do have torsion, the local exponential curve fits 
do not need to carry torsion in order to achieve good 
results in the local frame adaptation, see e.g. Fig. [ 8 ] 

The constraint of zero torsion forces us to split our 
exponential curve fit into a two-fold algorithm: 

Step 1 Estimate at g G SE( 3) the spatial velocity part 
c ( 1 )(^f) from the spatial structure tensor. 

Step 2 Move to a different location g new G SE( 3) 
where a horizontal exponential curve fit makes sense 
and then estimate the angular velocity c^ 2 ^ from the 
rotation part of the structure tensor over there. 

This forced splitting is a consequence of the next lemma. 


Lemma 4 Consider the class of exponential curves with 
nonzero spatial velocity ^ 0 such that their spatial 
projections do not have torsion. Within this class the 
constraint c 6 = 0 does not impose constraints on cur¬ 
vature if and only if the exponential curve is horizontal. 


Proof For a horizontal curve 7 g(t) we have x = 0 ^ 

c ( 1 ) -c (2) |«| 


c 1 = c 2 = 0 and indeed |r| = c || c ^i)n '^ 1 = c b \k\ = 0 
and we see that constraints c 6 = 0 and |r| = 0 reduce 
to only one constraint. The curvature magnitude stays 
constant along the exponential curve and the curvature 
vector at t = 0, recall Eq. (21), is in this case given by 


k( 0) = 1 -c 4 c 


c 5 c 3 

4^3 


is both orthogonal to vector = (c 1 , c 2 , c 3 ) T and or¬ 
thogonal to (—c 2 ,c 1 , 0 ) T , and thereby constrained to a 
one dimensional subspace. □ 

From these observations we draw the following conclu¬ 
sion for our exponential curve fit algorithms. 

Conclusion: In order to allow for all possible cur¬ 
vatures in our torsion-free exponential curve fits we 
must relocate the exponential curve optimization at g G 
SE(3) in U : SE( 3) -A M to a position g new G SE( 3) 
where a horizontal exponential curve can be expected. 
Subsequently, we can use Lemma [3] to transport the 
horizontal and torsion-free curve through g new , back to 
a torsion-free exponential curve through g. 

This conclusion is the central idea behind our fol¬ 
lowing two-fold algorithm for exponential curve fits. 


Algorithm Two-fold Approach: 


The algorithm follows the subsequent steps: 

Step la: Initialization. Compute structure tensor 
S s,p (g) from input image!/ :R 3 xS 2 -A R + via Eq. (97). 
Step lb: Find the optimal spatial velocity: 


C m (g) = argmin |^ C Q ^ M^S s ’ p (g)M fJ ,2 
l|c (1) || = M -1 


( 100 ) 


for g = (x, R n ), which boils down to finding the eigen¬ 
vector with minimal eigenvalue of the 3x3 spatial sub¬ 


matrix of the structure tensor (89). 


Step 2a: Given (g) we aim for an auxiliary set of 
coefficients, where we also take into account rotational 
velocity. To achieve this in a stable way we move to a 
different location in the group: 


Qnew — (x, Rn n 




= RnC 


( 1 ) 


( 101 ) 


and apply the transport of Lemma [3] afterwards. At 
g ne w , we enforce horizontally, see Remark [13] below, 
and we consider the auxiliary optimization problem 


which can be any vector orthogonal to spatial velocity 
c^ 1 ) = (0,0,c 3 ) T . Now let us check whether the condi¬ 
tion is necessary. Suppose t i-A 7 ®(t) is not horizontal, 
and suppose it is torsion free with c 6 = 0. Then we have 
(Ac 4 + c 2 c 5 = 0 , as a result the initial curvature 


K (°) = ETi 


C new(gnew) = argmin {c T M /i 2 S s ’ p (^ new; )M / , 2 c}. 
c e m 6 , 


IMG = 1, 

c 1 =c 2 =c 6 =0 


( 102 ) 


Here zero deviation from horizontality (34) and zero 
torsion ( 22 ) is equivalent to the imposed constraint: 

X = 0 and \t\ = 0 c 1 = c 2 = c 6 = 0 . 
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Step 2b: The auxiliary coefficients c new (g new ) = 
( 0 ,0, c 3 (g new ), c A (g new ), c 5 (g new ), 0 ) T of a torsion-free, 
horizontal exponential curve fit 7 ®”®“ through g new . 
Now we apply transport (via Lemma [3| of this hori¬ 
zontal exponential curve fit towards the corresponding 
exponential curve through g: 

c} inai ( 5 )= R r° )c new (g new ). (103) 

\ n new / 


This gives the final, torsion-free, exponential curve fit 
t jg ^ 9 \t ) in SE(3 ), yielding the final output pro¬ 
jected curve fit 

t ^ (7 © (0,e z ) € K 3 x S 2 , (104) 


with g = (x, R n ), recall Eq. (74). 


Remark 13 In step 2a of our algorithm we jump to a 
new location g new = (x, H nriew ) with possibly different 
orientation n new such that spatial tangent vector 


3 


E C * ^l(x,R„ nera ) > 

i= 1 


points in the same direction as n new G 5 2 , recall Eq. (31), 
from which it follows that n new is indeed given by (101). 
If (4 1 ) = (c 1 , c 2 , c 3 ) T = a = (0,0,1) T then n new = n. 


Lemma 5 The preceding algorithm is well-defined on 
the quotient M 3 x S 2 = SE{3) / ({0} x 50(2)). 


Proof To show that the preceeding algorithm is well- 
defined on the quotient we need to show that the final 
result ( fl04| is independent on both the choice of of 
R n G 50(3) s.t. H n e z = n and the choice of R nrieu; £ 
50(3) s.t. R Unew e z = n new . 

First, we show independence on the choice of R n . 

We apply Lemma [ 2 ] and thereby it remains to be shown 
that condition c } inal (gh a ) = Z^c } inal (g) is satisfied. 
This follows directly from Eq. ( |103[ ) if as long as n new 
found in Step 2 a is independent of the choice of R n . 
This property indeed follows from c^ 1 ) (gh a ) = R<T (g) 
which can be proven analogously to (J99| . Then we have 


j{gh a ) — RnR< 
= R n R< 


di) 


Ri „c W(g) 


C 

T 

'e 2 ,a AL e 2 ,a 


(gha) 


n„ 


,(g)- 


(105) 


So we conclude that (104) is indeed independent on the 
choice of R n . 


Finally, Eq. (104) is independent of the choice of 


Hn new - This follows from c new (gh a ) = Z a c new (g) in 


Step 2a. Then c^ inal in Eq. (103) is independent of the 
choice of R nri , ett; because Z J in c new Z^c new is can¬ 
celed by R, 


1 ^ Rn ew Re 


in Eq. (103). 


□ 



Fig. 8 : Volume rendering of a 3D test-image. The 
real part of the orientation score (cf. Section 7.1) 
provides us a density U on M 3 x 5 2 . Left: spatial 
parts of exponential curves (in black) aligned with 
spatial generator *4 3 |( xR r q, where n rnax (x.) = 


n max ( x ) 

argmax|[/(x, n)|. Right: spatial parts of our exponen- 

nG-S' 2 

tial curve fits Eq. (104) computed via the algorithm in 
Section 6.3.3[ which better follow the curvilinear struc¬ 
tures. 


In Fig. [ 8 ] we provide an example of spatially projected 
exponential curve fits in SE( 3) via the twofold ap¬ 
proach. Here we see that the resulting gauge frames 
better follow the curvilinear structures of the data (in 
comparison to the normal left-invariant frame). 


6.4 Exponential Curve Fits in SE( 3) of the 2 nd Order 

In this section we will generalize Theorem [2] to the case 
d = 3, where again we include the restriction to torsion- 
free exponential curves. 

6.4-1 The Hessian on SE(3) 

For second order curve fits we consider the following 
optimization problem: 


c *(#) — argmin 

ceM 6 ,||c|| -l,c 6 =0 




t =0 


(106) 


with V = G S *C/. Before solving this optimization prob¬ 
lem in Theorem[ 6 ]we first define the 6 x 6 non-symmetric 
Hessian matrix by 

(H B (U))(g) = [A j A i (V)](g), with V = G s * U (107) 

and where i = 1,...,6 denotes the row index, and 
j = 1, —, 6 denotes the column index. Again we write 
H S :=H S (V). 

Theorem 6 (Second Order Fit via Symmetric Sum 
Hessian) Let g G SE( 3) be such that the symmetrized 
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Hessian matrix \M~ 1 (H s (g)-\-(H s (g)) T )M ~ 1 has eigen¬ 
values with the same sign. Then the normalized eigen¬ 
vector M p c*(g) with smallest absolute non-zero eigen¬ 
value of the symmetrized Hessian matrix provides the 
solution c*(g) of optimization problem (106 1 ). 

Proof Similar to the proof of Theorem [5] (only now 
with summations from 1 to 5). Again we include our 
additional constraint c 6 = 0 by taking the smallest non¬ 
zero eigenvalue. □ 


Remark 14 The restriction to g E SE( 3) such that the 
eigenvalues of the symmetrized Hessian carry the same 
sign is necessary for a unique solution of the optimiza¬ 
tion. Note that in case of our first order approach via 
the positive definite structure tensor, no such cumber¬ 
some constraints arise. In case g E SE( 3) is such that 
the eigenvalues of the symmetrized Hessian have differ¬ 
ent sign there are 2 options: 


1. Move towards a neighboring point where the Hes¬ 
sian eigenvalues have the same sign and apply trans¬ 
port (Lemma [3| Fig. [7]) of the exponential curve fit 
at the neighboring point. 

2. Take c *(g) still as the eigenvector with smallest ab¬ 
solute eigenvalue (representing minimal absolute prin¬ 
cipal curvature), though this no longer solves (106). 


6.4-2 Torsion-free Exponential Curve Fits of the 2nd 
Order via a Two-Fold Algorithm 


In order to obtain torsion-free exponential curve fits 
of the second order via our two-fold algorithm, we fol¬ 
low the same algorithm as in Subsection 6.3.3[ but now 
with the Hessian field H s (107) instead of the structure 
tensor field. 

Step la: Initialization. Compute Hessian 
H s (g) from input image U : M 3 x S' 2 —x M + via Eq. (107). 

Step lb: Find the optimal spatial velocity by (100) 
where we replace M M 2 S s,p (g)M M 2 by H s (g). 

Step 2a: We again fit a horizontal curve at g new 
given by (101). The procedure is done via ( |102| ) where 
we again replace M /i 2 S s,p> (^)M /x 2 by H s (g). 

Step 2b: Remains unchanged. We again apply 
Eq. (|Tofl) and Eq. (fl04). 

There are some serious computational technicalities 
in the efficient computation of the entries of the Hessian 
for discrete input data, but this is outside the scope of 
this article and will be pursued in future work. 


Remark 15 In Appendix[C]we propose another two-fold 
second-order exponential curve fit method. Here one 
solves a variational problem for exponential curve fits 
where exponentials are factorized over respectively spa¬ 
tial and angular part. Empirically, this approach per¬ 
forms good (see e.g. Fig.[9|. 



Fig. 9: In black the spatially projected part of expo¬ 
nential curve fits t i-A Jg(t) of the second kind (fitted 
to the real part of the 3D invertible orientation score, 
for details see Fig. 10) of the 3D image visualized via 
volume rendering. Left: Output of the 2-fold approach 
outlined in Subsection |6.4.2| Right: Output of the 2- 
fold approach outlined in Appendix [c] with s p = \, 
s 0 = i(0.4) 2 , /I = 10. 


7 Image Analysis Applications 

In this section we present examples of applications where 
the use of gauge frame in SE(d ) obtained via expo¬ 
nential curve fits is used for defining data-adaptive left 
invariant operators. Before presenting the applications, 
we start by briefly summarizing the invertible orienta¬ 
tion score theory in Sec. |7.1| 

In case d = 2 the application presented is the en¬ 
hancing of the vascular tree structure in 2D retinal im¬ 
ages via differential invariants based on gauge frames. 
This is achieved by extending the classical Frangi ves- 
selness filter [37 to distributions U on SE( 2). Gauge 
frames in SE( 2) can also be used in non-linear multiple- 
scale crossing preserving diffusions as demonstrated in 
[66], but we will not discuss this application in this pa¬ 
per. 

In case d = 3 the envisioned applications include 
blood vessel detection in 3D MR-angiography, e.g. the 
detection of the Adamkiewicz vessel, relevant for surgery 
planning. Also in extensions towards fiber-enhancement 
of diffusion-weighted MRI [311130] the non-linear diffu¬ 
sions are of interest. Some preliminary practical results 
have been conducted on such 3D-datasets [44 1 124 1 122] . 
but here we shall restrict ourselves to very basic artifi¬ 
cial 3D-datasets to show a proof of concept, and leave 
these three applications for future work. 


7.1 Invertible Orientation Scores 

In the image analysis applications discussed in this sec¬ 
tion our function U \ M, d x i S d ~ l R is given by the 
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Fig. 10: Visualization of one iso-level of the real part 
of an invertible orientation score of a 3D-image created 
via the cake wavelet in Fig. 


11 


real part of an invertible orientation score: 

U (x, n) = Re{W^,/(x, R n )}, 

where R n is any rotation mapping reference axis a onto 
n G 5 d_1 , where / G L 2 (M d ) denotes a input image, and 
where ^ is a so-called ’cake-wavelet’ and with 

WV/(x, R n ) = f '0(Rn 1 (y - x))/(y) dy. (108) 

For d > 2 we restrict ourselves to wavelets satisfying 


^(Ra,a x ) = ^( x ) 5 f° r a ii X G R d (109) 

and for all rotations R aQ , €= Stab (a) (for d = 3 this 
means for all rotations about axis a, Eq. pi). As a re¬ 
sult U is well-defined on the left cosets x S d_1 = 
SE(d)/({ 0} x SO(d— 1)) as the choice of R n G SO(d ) 
mapping a onto n is irrelevant. See Fig. [lO] for an ex¬ 
ample of a 3D orientation score. 

If we restrict to disk-limited images, exact recon¬ 
struction is performed via the adjoint: 


/ = = t - 1 


& ^ —d - 

(2tt)2 


f -F[V^/(-, R)] ( W )^(R- 1 a>)d nso { d) (R) 

SO(d) 


( 110 ) 



Fig. 11: Visualization of cake-wavelets in 2D (top) and 
3D (bottom). In 2D we fill up the ‘pie’ of frequencies 
with overlapping “cake pieces”, and application of an 
inverse DFT (see m provides wavelets whose real and 
imaginary parts are respectively line and edge detec¬ 
tors. In 3D we include anti-symmetrization and the 
Funk transform [23 on I^-S 2 ) to obtain the same, see 
M- The idea is to redistribute spherical data from ori¬ 
entations towards circles laying in planes orthogonal to 
those orientations. Here, we want the real part of our 
wavelets to be line detectors (and not plate detectors) 
in the spatial domain. In the figure one positive iso¬ 
level is depicted in orange and one negative iso-level is 
depicted in blue. 


if ip is an admissible wavelet. The condition for admissi¬ 
bility of wavelets ip are given in [25] . In this article, the 
wavelets ip are given either by the 2D ‘cake-wavelet s’ 
used in mM or by their recent 3D-equivalents given 
in [44] . Detailed formulas and recipes to construct such 
wavelets efficiently can be found in [44 and in order to 
provide the global intuitive picture they are depicted in 

Fig-P1 

In the subsequent sections we consider two types 
of operators acting on the invertible orientation scores 
(recall d> in the commuting diagram of Fig.[3|: 

1. for d = 2, differential invariants on orientation scores 
based on gauge frames {£>i,£> 2 ,# 3 }- 

2. for d = 2,3, non-linear adaptive diffusions steered 
along the gauge frames, i.e. 


W (x, n, t) = W (x, R n , t) = * t (C0(x,R n ), (111) 
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where W(g,t ), with t > 0, is the solution of: 

{ rid 

a -^(g,t)=EDu(B i )\W(g,t), (n2) 

W(g,0) = U(g), 

where the gauge frame is induced by an exponential 
curve fit to data U at location g E SE(d). 

7.2 Experiments in SE( 2) 


Similarly to the vesselness filter , we need a mech¬ 
anism to robustly deal with vessels of different width. 
This is why for this application we extend the (all¬ 
scale) orientation scores to multiple-scale invertible ori¬ 
entation scores. Such multiple-scale orientation scores 
[66] coincide with wavelet transforms on the similitude 
group SIM(2) = M 2 x SO( 2) x M + , where one uses a B- 
spline |71I55] basis decomposition along the log-radial 
axis in the Fourier domain. In our experiments we used 
N = 4,12 or 20 orientation layers and a decomposition 
centered around M = 4 discrete scales ai given by 


We consider the application of enhancing and detecting 
the vascular tree structure in retinal images. Such im¬ 
age processing task is highly relevant as the retinal vas¬ 
culature provides non-invasive observation of the vas¬ 
cular system. A variety of diseases such as glaucoma, 
age-related macular degeneration, diabetes, hyperten¬ 
sion, arteriosclerosis or Alzheimer’s affect the vascula¬ 
ture and may cause functional or geometric changes 
[43] . Automated quantification of these defects promises 
massive screenings for vascular-related diseases on the 
basis of fast and inexpensive retinal photography. To 
automatically assess the state of the retinal vascular 
tree, vessel segmentation are needed. Because retinal 
images usually suffer from low contrast on small scales, 
the vasculature in the images needs to be enhanced 
prior to the segmentation. One well-established approach 
is the Frangi vesselness filter [37] which is used in robust 
retinal vessel segmentation methods |15il53| . However, 
a drawback of the Frangi filter is that it can not han¬ 
dle crossings or bifurcations that make up an important 
part of the vascular network. This is precisely where the 
orientation score framework and the presented locally 
adaptive frame theory comes into play. 

The SE(2 )-vesselness filter, extending Frangi ves¬ 
selness [37 to SE( 2) (cf. [42]) and based on the locally 
adapted frame {Si, £> 2 ? # 3 } is given by the following left 
invariant operator: 


$(U) = ’ e 
' 0 


(*- 


2<7 2 ^ if 0 > 0 , 

if a < 0 . ’ 

with anisotropy measure: 93 = ^ 7 

structureness: 6 = (B^U) 2 + (b\\J EB^IJ) 2 
convexity: 0 = B 2 U + B 2 U , 


(113) 


with <ji = \ and 02 = 0.2||S 2 2 ^ EB 2 U\\oq. Here the de¬ 
composition of the vesselness in structureness, anisotropy 
and convexity follows the same general principles of the 
vesselness. As in vessels are line like structures we use 
the exponential curve fits of 2nd order obtained via the 
symmetric product of the Hessian (i.e. solving the op¬ 
timization problem in Thm. [3|. 


ai = CLm^e 1 ( M 1 lo g(a maa; / a min ) ? (114) 

l = 0,..., M — 1 where a ma x is inverse proportional 
to the Nyquist-frequency p n and a m in close to the in¬ 
ner scale [34. induced by sampling (see [66] for details). 
Then, the multiple-scale orientation score is given by 
the following wavelet transform Wpf : SIM{2) -A C: 

W</,/(x,0,a) =* f ^(a _1 Rg 1 (y — x))/(y) dy, (115) 

R d 


and we again set U := Re{W^/}. Finally we define the 
total integrated multiple scale SIM (2)-vesselness by: 


(^ 57M (2)(t/))( x ) := 

M—l N 

Moo 1 E mA EWnADKMi), 

i=0 j =1 


(116) 


where SE(2 )-vesselness operator d> is given by Eq. (113), 
and where and denote maxima w.r.t. sup-norm 
|| • ||oo taken over the subsequent terms. 

Note that another option for constructing a SIM(2)~ 
vesselness is to use the non-adaptive left-invariant frame 
{4i, 4-2, * 4 . 3 } instead of the gauge frame. This non-adaptive 
SE(2 )-vesselness operator is obtained by simply replac¬ 
ing the Bi operators by the Ai operators in Eq. (113) 
accordingly. 

The aim of the experiments presented in this section 
is to show the following advantages: 


Advantage 1: The improvement of considering the multiple- 
scale vesselness filter via gauge frames in SE( 2), 
compared to multiple-scale vesselness [37] acting di¬ 
rectly on images. 

Advantage 2: Further improvement when using the 
gauge frames instead of using the left-invariant vec¬ 
tor fields in SE(2 )-vesselness (113). 


In the following experiment, we test these 3 tech¬ 
niques (Frangi vesselness [37] . SIM (2)^-vesselness via 
the non-adaptive left invariant frame, and the newly 
proposed SIM(2 )-vesselness via gauge frames) on the 








22 


Duits, Janssen, Hannink, Sanguinetti 



Fig. 12: From left to right: Retinal image / (from HRF- 
database), multi-scale vesselness filtering results for the 
multi-scale Frangi vesselness filter on M 2 , our SIM {2/ 
vesselness via invertible multi-scale orientation score 
based on left-invariant frame {Ai, A 2 , ^. 3 }, and based 
on adaptive frame {B 1 , 62 ,£> 3 }. 




publically availably High Resolution Fundus (HRF)- 
dataset [49], containing manually segmented vascular 
trees by medical experts. The HRF-dataset consists of 
wide-held fundus photographs for a healthy, diabetic 
retinopathy and a glaucoma group (15 images each). A 
comparison of the 3 vesselness filters on a small patch is 
depicted Fig.[l 2 | Here, we see that our method performs 
better both at crossing and non-crossing structures. 

To perform a quantitative comparisson, we devised 
a simple segmentation algorithm to turn a vesselness 
filtered image V(/) into a segmentation. First an adap¬ 
tive thresholding is applied, yielding a binary image 

f B m 6 >([V(/) — G 7 * V(/)] — h), (117) 


where 0 is the unit step function, G 1 is a Gaussian 
of scale 7 = \cj 2 1 and h is a threshold param¬ 
eter. In a second step, the connected morphological 
components in fs are subject to size and elongation 
constraints. Components counting less than r pixels or 
showing elongations below a threshold v are removed. 
Parameters 7 , r and v are fixed at 100 px, 500 px and 
0.85 respectively. The vesselness map V(/) : M 2 —>> M is 
one of the 3 methods considered. 

The segmentation algorithm described above is eval¬ 
uated on the HRF dataset. Average sensitivity and ac¬ 
curacy over the whole dataset are shown in Fig. [I3| as 
a function of the threshold value h. It can be observed 
that our method performs considerably better than the 
one based on the multi-scale Frangi filter. The segmen¬ 


tation results obtained with SIM ( 2 )-vesselness (116) 
based on gauge frames are more stable w.r.t variations 
in the threshold h and the performance on the small 
vasculature has improved as measured via the sensitiv¬ 
ity. Average sensitivity and accuracy at a threshold of 
h = 0.05 compare well with other segmentation meth¬ 
ods evaluated on the HRF dataset for the healthy cases 
(see m Tab. 5] and [42]). On the diabetic retinopa- 


3 cf. http://www5.cs.fau.de/research/data/fundus-images/ 


Fig. 13: Left: Comparison of multiple scale Frangi ves¬ 
selness and SIM ( 2 )-vesselness via gauge frames. Aver¬ 
age accuracy and sensitivity on the HRF dataset over 
threshold values h. Shaded regions correspond to ±1 a. 
Right: comparing of SIM ( 2 )-vesselness with and with¬ 
out including the gauge frame (i.e. using {Ai,A 2 ,As} 
in Eq. (fTl3|)). 


thy and glaucoma group, our method even outperforms 
existing segmentation methods. 

Finally, regarding the second advantage we refer to 
where the SIM ( 2 )-vesselness-filtering via the 
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Fig. 

locally adaptive frame produces a visually much more 
appealing soft-segmentation of the blood vessels than 
SIM (2)-vesselness filtering via the non-adaptive frame. 
It therefore also produces a more accurate segmenta¬ 
tion as can be deducted from the comparison in Fig .[13] 
For comparison, the multiscale Frangi vesselness filter 
is also computed via summation over single scale results 
and max-normalized. Generally, we conclude from the 
experiments that the locally adaptive frame approach 
better reduces background noise, showing much less 
false positives in the final segmentation results. This 
can be seen from the typical segmentation results on 
relatively challenging patches in Fig. [15] 


7.3 Experiments in SE(3) 


We now show first results of the extension of coherence 
enhancing diffusion via invertible orientation scores (CE 
DOS [36] ) of 2 D images to the 3D setting. Again, data 
is processed according to Fig. [3] Fir st, we construct an 
orientation score according to ( |108| ), using the 3D cake 
wavelets (Fig.[TT]). For determining the gauge frame we 
use the first order structure tensor method in combi¬ 
nation with Eq. (118) in Appendix [A] In CEDOS we 
have df> = <f> t ^ as defined in ( 111 ) and ( | 112 ), which is a 
diffusion along the gauge frame. 
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Fig. 14: Center: original image from HRF-dataset (healthy subject nr. 5). Rows: the soft-segmentation (left) and 
the corresponding performance maps (right) based on the hard segmentation (117). In green: true positives, in 
blue true negatives, in red false positives, compared to manual segmentation by expert. 1 st row: SIM( 2 )-vesselness 
(116) based on non-adaptive frame {Ai, A 2 , A 3 }. 2 nd row: SIM( 2 )-vesselness (116) based on the gauge frame. 


Image-part Manual exp. Non-gauge Gauge 



Fig. 15: Two challenging patches (one close to the op¬ 
tic disk and one far away from the optic disk processed 
with the same parameters). The gauge frame approach 
typically reduces false positives (red) on the small ves¬ 
sels, and increases false negatives (blue) at the larger 
vessels. The top patch shows a missing hole at the top 
in the otherwise reasonable segmentation by the expert. 


The diffusion in CEDOS can enhance elongated struc¬ 
tures in 3D data while preserving the crossings as can 
be seen in the two examples in Fig.[l 6 j In these exper¬ 
iments as well as in the example used in Fig. [ 8 ] [9] and 
[To] we used the following 3D cake-wavelet parameters 
for constructing the 3D-invertible orientation scores: 
N 0 = 42,5 0 = 0.7, k = 2,7V = 20,7 = 0.85, L = 16 
evaluated on a grid of 21 x 21 x 21 pixels, for details see 


[44] . The settings for tangent vector estimation using 
the structure tensor are s p = |(1.5) 2 ,5o = 0 , p Q = 
|(0.8 ) 2 and p = 0.5. We used p p = ^( 2) 2 for the 
first dataset (Fig. [l 6 | top), and p p = ^(3.5 ) 2 for the 
second dataset (Fig. |16| bottom). For the diffusion we 
used t 2.5, .Du = D 22 0 . 01 , D33 = 1 ,D 44 = 

D 55 = D 66 = 0.04, where the diffusion matrix is given 
w.r.t. gauge frame {61,62, # 3 , 64, 65, Bo}, and normal¬ 
ized frame {p~ x A\, p~ x A2, p~ x A3, A4, A3, Aq} . 


The advantages of including the gauge frames w.r.t. 
the non adaptive frame can be better appreciated in 
Fig-[l3 Here, we borrow from the neuroimaging com¬ 
munity the glyph visualization, a standard technique for 
displaying distributions U : M 3 x S 2 -A M + . In such vi¬ 
sualizations every voxel contains a spherical surface plot 
(a glyph) in which the radial component is proportional 
to the output-value of the distribution at that orienta¬ 
tion, and the colors indicate the orientations. One can 
observe that diffusion along the gauge frames include 
better adaptation for curvature. This is mainly due to 
the angular part in the S 3 -direction, cf. Fig. [18] which 
includes curvature, in contrast to 4 . 3 -direction. The an¬ 
gular part in S 3 causes some additional angular blurring 
leading to more isotropic glyphs. 
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(c) Curve fits 



(d) Slice+Noise (e) Gauge (f) No Gauge 





(j) Slice+Noise 


(k) Gauge (1) No Gauge 


Fig. 16: Results of CEDOS with and without the use 
of gauge frames, on 3D artificial datasets containing 
highly curved structures. Gauge frames are obtained, 
Appendix [Aj via 1st order exponential curve fits 


using the two-fold algorithm of Subsection 6.3.3 


8 Conclusion 

Locally adaptive frames (‘gauge frames’) on images based 
on the structure tensor or Hessian of the images are ill- 
posed at the vicinity of complex structures. Therefore 
we create locally adaptive frames on distributions on 
SE(d ), d = 2,3 that extend the image domain (with 
positions and orientations). This gives rise to a whole 
family of local frames per position, enabling us to deal 
with crossings and bifurcations. In order to generalize 
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(c) Enhanced using frame (d) Enhanced using frame 
{Hi,...,5 6 } {Ai,...,As} 


Fig. 17: Glyph visualization (see text) of the absolute 
value of the diffused orientation scores with and with¬ 


out the use of gauge frames in the the artificial dataset 


depicted in Fig. 16 top. 


gauge frames in the image domain to gauge frames in 
SE(d ), we have shown that exponential curve fits gives 
rise to suitable gauge frames. We distinguished between 
exponential curve fits of the 1st order and of the 2nd 
order: 

1. Along the 1st order exponential curve fits, the 1st 
order variation of the data (on SE(d)) along the ex¬ 
ponential curve is locally minimal. The Euler-Lagrange 
equations are solved by finding the eigenvector of 
the structure tensor of the data, with smallest eigen¬ 
value. 

2. Along the 2nd order exponential curve fits, a 2nd or¬ 
der variation of the data (on SE(d)) along the expo¬ 
nential curve is locally minimal. The Euler-Lagrange 
equations are solved by finding the eigenvector of 
the Hessian of the data, with smallest eigenvalue. 

In SE( 2), the 1st order approach is new while the 2nd 
order approach formalizes previous results. In SE( 3), 
these two approaches are presented for the first time. 
Here, it is necessary to include a restriction to torsion- 
free exponential curve fits in order to be both compati¬ 
ble with the null-space of the structure/Hessian tensors 
and the quotient structure of M 3 xi S 2 . We have pre¬ 
sented an effective two-fold algorithm to compute such 
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torsion-free exponential curve fits. Experiments on ar¬ 
tificial datasets show that even if the elongated struc¬ 
tures have torsion, the gauge frame is well-adapted to 
the local structure of the data. 

Finally, we considered the application of a differ¬ 
ential invariant for enhancing retinal images. Experi¬ 
ments show clear advantages over the classical vessel- 
ness filter [37] , Furthermore, we also show clear ad¬ 
vantages of including the gauge frame over the stan¬ 
dard left-invariant frame in SE( 2). Regarding 3D im¬ 
age applications, we managed to construct and imple¬ 
ment crossing-preserving coherence enhancing diffusion 
via invertible orientation scores (CEDOS), for the first 
time. However, it has only been tested on artificial datasets. 
Therefore, in future work we will study the use of lo¬ 
cally adaptive frames in real 3D medical imaging appli¬ 
cations, e.g. in 3D MR angiography [45]. Furthermore, 
in future work we will apply the theory of this work and 
focus on the explicit algorithms, where we plan to re¬ 
lease Mathematica -implementations of locally adaptive 
frames in SE( 3). 
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A Construction of the Locally Adaptive Frame 
from an Exponential Curve Fit 

n d 

t ^2 c l Ai 

Let 7 g(t) = ge i=1 be an exponential curve through g 
that fits data U : SE(d) —» R at g G SE(d) in Lie group 
SE(d) of dimension n d — d(d + l)/2. In Section E(d = 2), 
and in Section[ 6 ](<i m 3), we provide theory and algorithms to 
derive such curves. In this section we assume 7 ^(-) is given. 

Recall from that the (physical) velocity at time t of 

n d 

the exponential curve 7 ^ equals ( 7 g)'(t) = c% Ai\ g= ~ c( ^ t y 

i= 1 9 

Recall that the spatial and respectively rotational compo¬ 
nents of the velocity are stored in the vectors 

c(!) = (c\...,c d ) T eM d , 

= (c d+1 ,..., c nd ) T 6 R 1 ^ 


Let us write c = 


c (i) 

cJ 2 ) 


with n d = d + r d . 


Akin to the case d = 2 discussed in the introduction we 
define the Gauge frame via B \= (R C ) T M“ 1 A, but now with 


B « (Bi,..., B nd ) T , An d ) T , 

M Id, 0 


“ ' 0 I, 


and R c = R 2 R 1 G SO{n d ). 


(118) 


For explicit formulae of the left-invariant vector fields in the 
d-dimensional case we refer to [32] . 

Now Ri is the counter-clockwise rotation that rotates the 

spatial reference axis reca ^ our conven Ii° n 

Ml|c (1) ||a" 
c ( 2 ) 

two vectors. Rotation R 2 is the counter-clockwise rotation 


onto 


strictly within the 2D-plane spanned by these 


that rotates 


Ml|c (1) ||a 


c (2) J onto ^ c (2) 

2D-plane spanned by these two vectors. As a result one has 


gc 


(i) 


strictly within the 


c = M- 1 R< 

H' 


R, /pHc^Ha 
^ 1 c ( 2 > 

a 
0 




yc 


(i) 


A2) 


(119) 


In particular we have that the preferred spatial direction 


• A is mapped onto 


0I -B = c-A. 


The next theorem shows us that our choice of assigning 
an entire gauge frame to a single exponential curve fit is the 
right one for our applications. 

Theorem 7 (construction of the gauge frame) Let c(g) de¬ 
note the local tangent components of exponential curve fit 
t 1 —y t ) at g = (x, R) G SE(d) in the data given by 

U(x,R) = U(x,Ra). Consider the mapping of the frame of 
left-invariant vector fields A\ g to the locally adaptive frame: 


B\g~(R <9) ) T M- 1 A\ g , (120) 

with R c = R 2 R 1 G SO{n d ), with s ubseq uent counter-clockwise 
planar rotations R \, R 2 given by [H9 1 ). Then the mapping 
A g B g has the following properties: 

— The main spatial tangent direction (L^)*^^ • A\ e is 

mapped to exponential curve fit direction c T (g) ■ A\ . 

— Spatial left-invariant vector fields that are 0 ^-orthogonal 
to this main spatial direction stay in the spatial part of 
the tangent space T g (SE(d)) under rotation R c and they 
are invariant up to normalization under the action \120\ 
if and only if the exponential curve fit is horizontal. 

Proof Regarding the first property we note that 

<£.).(;)• 4.-(;)• 4, 

as left-invariant vector fields are obtained by push-forward of 
the left multiplication. Furthermore, by Eq. ( |120| ) and Eq. (U9} 
we have 

(o)-B = M- 1 R c (^)-^ = c.A 
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Fig. 18: Visualization of the mapping of left-invariant 
frame {Ai, ..., ^4.6} | ^ onto locally adaptive spatial frame 
{23i,..., & 6 }\ g and 7^(*) a non-horizontal and torsion-free 
exponential curve passing trough g = e = (0, 1). The top row 
indicates the spatial part in R 3 , whereas the bottom row in¬ 
dicates the angular part in S 2 . The top black curve is the 
spatial projection of 7 g(-), and the bottom black curve is the 
angular projection of the exponential curve. After application 
of the exponential curve is horizontal w.r.t. the frame 
{Ai}, subsequently leaves spatial generators orthogonal 
to a horizontal curve invariant, so that A\ = B\ and A2 = B2 
are strictly spatial, as is in accordance with Theorem [7] The 
angular part of B3 is shown via the curvature at the blue ar¬ 
row. The spatial part of B 4 and £>5 is depicted in the center 
of the ball. 

Regarding the second property, we note that if b • a = 0 =>- 

and 7 g is horizontal iff || = a in which case the planar ro¬ 
tation R 2 reduces to the identity and R 2 R 1 = (o) 

and only spatial normalization by /r _1 is applied. □ 

Remark 16 For d = 2 and a = (1,0) T the above theorem 
can be observed in Fig.[5] where main spatial direction Ai = 
cos 6 d x + sin 6 d y is mapped onto B\ = c • A and where A2 
is mapped onto B 2 = sinyAi TcosyA^)- 

Remark 17 For d = 3 and a — (0,0,1) T the above theorem 
can be observed in Fig. |18| where main spatial direction A3 = 
n • V^s is mapped onto B3 = c • A, and where Ai and A2 
are mapped to the strictly spatial generators B\ and B2 • For 
further details see m- 


B The Geometry of Neighboring Exponential 
Curves 

In this appendix we provide some differential geometry un¬ 
derlying the family of neighboring exponential curves. 


First we prove Lemma [3] on the construction of the fam¬ 
ily { 7 ? g } of neighboring exponential curves in SE( 3), recall 
Fig. [7] and then we provide an alternative coordinate free 
definition of 7 ^ g in addition to our Definition |5j . 

For the proof of Lemma [3] we will just show equalities 
( | 86 | ) as from this equality it directly follows by differentiation 
w.r.t. t that the exponential curves 7 ^ g {-) = (x^(-), R^(-)) 
and 7 g(-) = (xg(-), R^(-)) have the same spatial and angular 
velocity. For the spatial velocities it is obvious, for the angular 
velocities, we note that rotational velocity matrices flh and 
fig are indeed equal: 

R hit) = R^(t)R — 1 R / => 
n h := ^R^W | t=0 (R-fi(o )) -1 

= ii n 9(t) | t= 0 (R B (o))- 1 -« B , 
where we note that R g (t) = e l ^>R 9 ( 0 ) = e*^«R, and 

R h (t) = e'^R^O) = e^R'. 

Regarding the remaining derivation of ( | 86 | ), we note that 
it is equivalent to 

7 n, g (t) = h (0, (R') T R) (g- 1 ^)) ( 0 , (R') t R)-^ (121) 

by group product 0 - So we focus on the derivation of this 
identity. Let us set Q := (R') T R. Now relying on the ma¬ 
trix representation 0 and matrix exponential we deduce the 
following identity for h = e = ( 0 , /): 

( Q °) c 

7 9 W = V° Q7 (t) = (0, Q)7e(t) (0,Q _1 ), (122) 

which holds for all Q G 50(3), in particular for Q = (R') T R. 
Consequently, we have 

7 h,g(t) = hX t C t,g{t) = h(0,Q)g~ 1 lg(t) (0, Q _1 ), 

from which the result follows. □ 

We conclude From Lemma [3] that our Definition [5] is indeed 
the right definition for our purposes, but as it is a definition 
expressed in left-invariant coordinates it also leaves the ques¬ 
tion what the underlying coordinate-free unitary map from 
T g (SE( 3)) to Th(SE(3)) actually is. Next we answer this 
question where we keep Eq. pit in mind. 

Definition 6 Let us define the unitary operator 
a h ’ B ■ Tg(SE( 3)) -»■ T h {SE{ 3)) by 

iX h ’ 8 ■= (L h ).R h -i g (L g -i)„ 

for each pair g = (x, R), h = (x', R') G SE(3). 

Remark 18 From \ 87 \ it follows that the unitary correspon¬ 
dence between T~ C ( t ) and T~c ( t ) is preserved for all t G R. 

Definition 7 The coordinate free definition of 7 ^ g is that 
it is the unique exponential curve passing through In at t = 0 
with 

(7h, s ) , (0) = iX h ’ 9 ((77(0)). 

Remark 19 The previous definition (Definition [ 5 ]) follows from 
the coordina te fre e definition (Definition]?]). This can be shown 
via identity (m} which can be rewritten as 

7h >9 (*) ~L h o conj(0, Q) o (L g -i)^(t)), (123) 

which indeed yields 

(L h o conj(0, Q) o (L r i))* = (L h )* o Ad(0, Q) o (L g - 1 )* 

= (^)*(? q)( l 9 -o* 

= U h ' 9 , 

again with Q = (R') T R, conj(g)h — ghg -1 , and Ad (g) = 
(conj(g))* is the adjoint representation m- 
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C Exponential Curve Fits on SE( 3) of the 2nd D The Hessian induced by the left Cartan 
Order via Factorization connection 


Instead of applying a 2nd order exponential curve fit ( |106| ) 
containing a single exponential one can factorize exponentials, 
and consider the following optimization: 


In this section we will provide a formal differential geometrical 
underpinning for our choice of Hessian-matrix 

H (17) = [Aj(AiU)\, (127) 


c *(g) = argmin 

cGM 6 ,||c|| m = 1,c 6 = 0 

I ^V( S e*( clAl + c2A 2 + c 3 A 3 ) e t( c 4 A 4 +C =A 5 )) 


(124) 


As shown in Theorem [8] the Euler-Lagrange equations are 
solved by spectral decomposition of the symmetric Hessian 
given by 


f AiAiV ... A!A 6 V\ 
W-.= M s (u)= : *, : 


yAiAeV • • • *46*46 E J 


(125) 


with V = G s * U . This Hessian differs from the consistent 
Hessian in Appendix [B] 


Theorem 8 (Second Order Fit via Factorization) Let 

g G SE(3) be such that Hessian matrix M~ 1 {H s {g))M~ 1 
has two eigenvalues with the same sign. Then the normalized 
eigenvector M^c*(g) with smallest eigenvalue provi des t he 
solution c* (g) of the following optimization problem \12J$ . 

Proof Define F ± := c^ 1 ) • G T e (SE( 3)) with := 
(A 1 ,A 2 ,A 3 ) t . Define F 2 := c< 2 ) • A< 2 ) G T e (SE( 3)) with 
A^ 2 ) : = (A4, A5, Aq) t . Define vector fields Ei\ g := (L g )*Fi, 
T 2 \ g := ( L g )*F 2 . Then Then 

*V(ge tF >e tF *) I = lirn V( g e^e^)-2 V( g ) + V( a e-^e-^) 

dt J \t=0 h-t- 0 h 

= V(g) + T 2 T 2 V{ g) + 2F 1 X 2 V{g) 
= {c(g)) T U s (g)c(g). 


This follows by direct computation and the formula 
V{qe hF “) = V(q) + hF k V{q) + + 0(h 3 ), 


applied for (q = ge hFl , k = 2) and (q = g, k = 1). 

Therefore we can express the optimization functional as 


£( c y— | A 1 +c 2 A 2 -hc 3 A 3 ) e t(c 4 A4+c 5 A 5 ) 

= |c T H s (a)c|, 


t = 0 


(126) 


where i denotes the row-index and j the column index on 
SE(d ), recall the case d = 2 in ( |62| ) and recall the case d = 3 in 
( |107| ). Recall from Theorem[2] Theorem[3]and Theorem[6]that 
this Hessian naturally appears via direct sums or products in 
our exponential curve fits of second order on SE(d). 

Furthermore we relate our exponential curve fit theory to 
the theory in [46] . where the same idea of 2nd order fits of 
auto-parallel curves to a given smooth function U : M —> R 
in a Riemannian manifold is visible in [461 Eq.3.3.50]. Here 
we stress that in the book of Jost [46t Eq.3.3.50] this is done 
in the very different context of the torsion-free Levi-Civita 
connections, instead of the left Cartan connection which does 
have non-vanishing torsion. 

Let us start with the coordinate free definition of the Hes¬ 
sian induced by a given a connection V* on the cotangent 
bundle. 


Definition 8 (coordinate free definition Hessian) On a Rie¬ 
mannian manifold (M, G) with connection V* on T*(M), the 
Hessian of smooth function U : M —> R is defined coordinate 
independently (|46) Def.3.3.5]) by V*dU. 

In coordinate-free form one has (cf. [461 Eq.3.3.50]) 


V*dU(X p , X p ) = A<W)) 


(128) 


for the auto-parallel (i.e. = 0) curve y(t) with tangent 

7 7 (0) = X p passing through c(0) = p at time zero. 

Remark 20 In many books on differential geometry V* is 
again denoted by V (we also did this in our previous works 
[29ll26| b In this appendix, however, we distinguish between 
the connection V on the tangent bundle and its adjoint con¬ 
nection V* on the co-tangent bundle T*(SE(d)). 

Let us recall that the structure constants of the Lie alge¬ 
bra are given by 


rid 

[Ai , Aj \ = AiAj — AjAi = c^Ak- (129) 

k=1 

As shown in previous work 29 the left Cartan connectior^ 
V on M — (SE(d), 0^), is the (metric compatible) connec¬ 
tion whose Christoffel symbols, expressed in the left-invariant 
moving (co)frame of reference, are equal to the structure con¬ 
stants of the Lie algebra: 

^•=4 = -^. e {- 1 , 0 , 1 }. 


with again boundary condition (p(c) = c T M 2 c = 1, from 
which the result follows via Euler-Lagrange V£ = A Vcp and 
left multiplication with M” 1 . □ 

This approach can again be decomposed in the two-fold 
approach. Effectively, this means that in Section |6.4.2| the up¬ 
per triangle of the Hessian H s is replaced by the lower trian¬ 
gle, whereas the lower triangle is maintained. This approach 
performs well in practice; see e.g. Fig. [9] where the results 
of the exponential curve fits of second order are similar to 
exponential curve fits of first order. 


More precisely, this means that if we compute the covariant 

Ud 

derivative of a vector field Y = y k Ak (i.e. a section in 

k =1 
rid . 

T(SE(d)) along the tangent y(t) m 1 % {^) -4^ly(t) °f some 

i = l 

smooth curve t 1 —>> y(t) in SE(d). This is done as follows 

n d ( ri d \ 

V T = E \ y k - E 4AV A h , (iso) 

k = 1 \ i,j = 1 J 


4 also known as minus Cartan connection 
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where we follow the notation in Jost’s book [461 p.108] and 
define y k (t) := -^y k {flit)). By duality this induces the fol¬ 
lowing (adjoint) covariant derivative of a covector field A (i.e. 
a section in T* (SE(d))): 

rid / ri d \ 

E A = E U + E ( 131 ) 

i = 1 \ k,j = 1 / 

with A i(t) = -^fXiflft)). Then by antisymmetry of the struc¬ 
ture constants it directly follows (see e.g. [26]) that the auto¬ 
parallel curves are the exponential curves: 

V -7 = 0 and 7 ( 0 ) = c and 7 ( 0 ) = g O 7 = 7 °. (132) 


Remark 21 Due to torsion of the left Cartan connection, the 
auto-parallel curves do not coincide with the geodesics w.r.t. 
metric tensor 0£. This is in contrast to the Levi-Cevita con¬ 
nection (see for example Jost’s book [46] ch:3.3, ch:4]) where 
auto-parallels are precisely the geodesics (see [461 ch:4.1]). 

Intuitively speaking this means that in the curved geome¬ 
try of the left Cartan connection on SE(d) (that is present in 
the domain of an orientation score, see Figure[3]) the ‘straight 
curves’ (i.e. the auto-parallel curves) do not coincide with the 
‘shortest curves’ (i.e. the Riemannian distance minimizers). 

The left Cartan connection is the consistent connection on 
SE(d) in the sense that auto-parallel curves are the expo¬ 
nential curves studied in this article. Therefore the consistent 
Hessian form on SE{d) is induced by the left Cartan connec¬ 
tion. Expressing it in the left-invariant frame yields 


V*d U(Ai,Aj) d = (V^d UflAfl 

(fToTT n d _ rid _ . 

£ (A i A r U+ E c^A k U)^\Aj) 

j' = 1 k = 1 

^ (AiAjU + E c%A k U) 

k = 1 

(Ai(AjU) + (AjAi — AiAj)U) 
_ - (Aj(AiU), _ 

(133) 

where i denotes the row-index and j the colu mn index. So 
we conclude from this computation that \127\ is the correct 
consistent Hessian on SE(d) for our purposes. 

Remark 22 The left Cartan connection has torsion and is not 
the same as the standard torsion-free Cartan-Schouten con¬ 
nection on Lie groups, which have also many applications in 
image analysis an statistics on Lie groups, cf. 1551551 ■ Recall 
that within the orientation score framework, right invariance 
is undesirable. 
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E Table of Notations 


Symbol 


Explanation 


Reference 


E.l Spaces and Input Data 


SE(d) 

The group of rotations and translations on R d 

Section 1 Section 2.1 and (4) 

M d xi S d 1 

Space of positions & orientations as a group quotient in SE(d) 


3 

), and Section 

6.1 


u 

Input data U : SE(d) — > R 


1 

), and Section 

2.1 


u 

Input data U : R d x S d ~ 1 — > R 


1 

), (3), and Section 6.1 

V 

Gaussian smoothed input data V = G s * U 

( 

2 ! 

61 , and (27) 


E.2 Tools from Differential Geometry 



Left-invariant vector field Ai restricted to g E SE(d) 

Section 

2.3 and (10), (9) 

Bil, 

Gauge vector field Bi restricted to g E SE(d) 

Section 

3| and (3( 

3) 

®mI 0 

Metric tensor restricted to g E SE(d) 

Section 

T5 

and ( 

23 


ll-IU 

//-norm on R nd , with rid = dim (SE(d)) = d ( d + 1 ^> 

Section 


and ( 

24 



Matrix := ( ^ 0 d J is used in definition of the /x-norm || • || M 

Section 

2.5 

and ( 

24 


d U(g) 

Derivative of U at g which is a covector in T*(SE(d)) 

( 12 ) 

vu{g) 

Gradient of U at g which is a vector in T g (SE(d)) 

Section 2.7 and L 

28) 


X 

Deviation from horizontality angle 

( 

34) 

c 

Left regular representation given by C g U{h ) = t/(g - 1 h) 

( 

7 

) 

71 

Right regular representation given by 7Z g U(h) = U(hg ) 

( 

7 

) 

L 

Left multiplication L g h — gh 

( 

8 

) 


E.3, part I: Exponential Curves and Exponential Curve Fits on SE(d) 


tS(-) 

Exponential curve starting from g with velocity c = (c^ 1 ),^ 2 )) 


16 

), and Section |2.4| 


Neighboring exponential curve starting at h E SE(d) with the same 
velocity as curve 7^ 

i 

51 

3.2 

) and Section |5.1| (|84|) and Section 

7 9 C ' (9, (-) 

Exponential curve fit to data U at g E SE(d) 


16 

), and Theorem 1|2 3|4|5|6 Fig. 3 

c * (d) 

Local tangent vector to exponential curve fit 7 ^* to data U 


54 

), (|63|), (64), ( 88 ), and (106) 

f 

Rotation in Th(SE(d)) arising in the construction of 7^ q 


52 

), and (85|) 

S S >P 

Structure tensor S s,p := S s,p (l7) of input data U 


91 

, (89), and (97) 

H s 

Gaussian Hessian H s := H s (t7) = H(V r ) of input data U 

( 

62), and (107) 


E.3, part II: Exponential Curve Fits on SE(3) with Projections in R 3 x 5 2 


Ra, (f> 

Counter-clockwise 3D rotation about axis a by angle cf> 

text below (31) 

Rn 

Any 3D rotation that maps a = (0, 0,1) T onto n E S 2 


31 

), and Theorem 5 

hex 

Element h a = (0,R a ,a) of the subgroup = {0} X 50(2) 


76 


O 

Symbol denoting action of SE(3) onto R 3 x S 2 


74 

), and Section 6.1 

Za 

Rotation matrix in 50(6) that arises in VC/ if g i—>■ gh a 

( 

80 

) 

M 

Null-space of the structure tensor S s,p 

( 

90 

) 

T'(y,n)(') 

Projected exponential curve fit to data U at (y, n) E R 3 x 5 2 


9j), and (98) 

gnew 

Location in SE{ 3) for horizontal exponential curve fit 


101 ) 


E.4 Applications 


/ 

Input greyscale image f : R d —>■ R 


108 

), and (110 

) 

MW 

Orientation score of greyscale image / via cakewavelet 'ip 


108 

), and Fig. 

3 10 11 


Nonlinear diffusion operator (diagonal diffusion in gauge frame) 


112 



Scale space representation of U at g E SE(d) and scale t > 0 


n2 


<7> 

Vesselness operator 


U3 



E.5 Appendix 


< 

Covariant derivative of vector field Y along 7 w.r.t. Left-Cartan 
connection V on T(SE(d)) 

(130 

), and Appendix jDj 

7 

Covariant derivative of covector field uj along 7 w.r.t. the adjoint 
Left-Cartan connection V* on T*(SE(d)) 

(131 

), and Appendix [d| 

V*dt7 

Coordinate-free definition of the Hessian 

(128 

i, < 

133), and Appendix D 
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